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ASYMPTOTIC THEORY OF TWO-DIMENSIONAL TRAILING- EDGE FLOWS* 

By R. E. Melnik and R. Chow 
Grumman Aerospace Corportation 

SUMMARY 

In this investigation, problems of laminar and turbulent viscous interaction near 
trailing e<^es of streamlined bodies are considered. The laminar study is based on the 
triple-deck formulation of Stewartson, Messiter, and Brown. This theory is developed 
from asymptotic expansions of the Navier- Stokes equations in the limit of large Reynolds 
numbers. The expansions describe the local solution near the trailing edge of cusped or 
nearly cusped airfoils at small angles of attack in compressible flow. A complicated 
inverse iterative procedure, involving finite- difference solutions of the triple-deck equa- 
tions coupled with asymptotic solutions of the boundary values, is used to accurately solve 
the viscous interaction problem. Results are given for the correction to the boundary- 
layer solution for drag of a finite flat plate at zero angle of attack. A solution is also 
presented for the viscous correction to the lift of an airfoil at incidence. A comparison of 
the present results with triple- deck solutions recently obtained by other investigators for 
the symmetric problem is presented. Also presented are some comparisons of the pres- 
ent solution with low Reynolds number (R s 200) solutions of the Navier- Stokes equations 
and with experimental data. These comparisons indicate that the asymptotic triple- deck 
theories are accurate over a surprisingly wide range of Reynolds numbers down to 
Reynolds numbers as low as 10 or less. 

In the second part of this investigation, the problem of turbulent interactions at air- 
foil trailing edges is considered. It is demonstrated that second-order boundary-layer 
theory fails at airfoil trailing edges and that the concept of the flow over an equivalent 
body formed from the displacement thickness is not appropriate for turbulent flows near 
trailing edges. A rational asymptotic theory is developed for treating turbulent interac- 
tions near trailing edges and is shown to lead to a multilayer structure of turbulent bound- 
ary layers. The flow over most of the boundary layer is described by a Lighthill model 
of inviscid rotational flow. The main features of the model are discussed and a sample 
solution for the skin friction is obtained and compared with the data of Schubauer and 
Klebanoff for a turbulent flow in a moderately large adverse pressure gradient. 


*This research was performed under NASA Langley Contract No. NAS 1-12426. 
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INTRODUCTION 


The problem of trailing- edge flows is of considerable importance in aerodynamics.' 
Most streamlined bodies end in a sharp trailing edge that is cusped or nearly cusped in ‘ 
order to provide a smooth transition of the flow into the wake. The flow near the trailing 
edge is important in establishing the lift and drag forces on the body. 

At high Reynolds numbers the solution of the Navier- Stokes equations can be 
expanded in descending powers (and logarithms) of the Reynolds number. The leading 
term is governed by inviscid flow equations over most of the domain and by boimdary- 
layer equations in a thin layer near the surface and in the wake. To lowest order, the lift 
is determined by solutions of the inviscid-flow equations, subject to the Kutta condition. 
Skin- friction drag is determined by solutions of the boundary-layer equationsi with the 
pressure distribution obtained from the inviscid solution. Within this approximation, 
form drag is computed from the surface pressures induced by the effect of the boundary 
layer on the external inviscid flow and is, therefore, a second- order effect in the theory. 

Although the inviscid and boundary- layer solutions provide the leading approxima- 
tion for the flow over streamlined bodies, higher order corrections are important in many 
problems. A major impediment in the determination of toe correction is due to the fact 
that the underlying asymptotic expansions are not uniformly valid at trailing edges. The 
nonuniformity is caused by the appearance of singularities in solutions of both toe laminar 
and turbulent boundary-layer equations at trailing edges. The nature of the singularity 
differs in toe laminar and turbulent cases, but in both cases, toe major effect is the pro- 
duction of a displacement thickness that is singular at the trailing edge. This in turn 
leads to singularities in toe induced pressures at the trailing edge. As a result, toe 
second- order Kutta condition cannot be satisfied and the viscous correction to lift cannot 
be determined. . In addition, corrections to the bovuidary- layer solutions for skin friction 
and form drag are not correctly determined by standard second-order boundary-layer 
theory. Most existing engineering methods for predicting viscous effects on lift are based 
on iterative solutions of the second- order boundary- layer equations. These methods expe- 
rience difficulties at trailing edges which are circumvented by an ad hoc smoothing of toe 
displacement surface determined from toe solution of toe boundary-layer equations. 

In spite of its importance and the continuing interest of many investigators, it is only 
recently that a comparatively complete theory of trailing-edge flows has been developed 
and this only for laminar flows. The recent advances in laminar trailing-edge problems 
are based on the triple- deck formulation of Stewartson (ref. 1) and Messiter (ref. 2) devel- 
oped originally for a symmetric flat plate at zero angle of attack. . These theories. were 
then extended the following year by Brown and Stewartson (ref. 3) to toe lifting flat plate. 
The triple-deck theories are applicable to general airfoils with a cusped or nearly. cusped 
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trailing edge in compressible flow. For Mach numbers not near one, compressibility 
enters, into the theory only through a scaling of dependent and independent variables. 
However, even though a number of properties of the solution were determined and the 
sin^ar behavior was explained, accurate numerical solutions were not obtained in these 
early works. 

The first part of this paper deals with the laminar- viscous interaction near airfoil 
trailing edges in the limit of large Reynolds numbers. The approach used is to develop 
the appropriate numerical procedures to solve the boundary-value problem formulated by 
Stewartspn, Messiter, and Brown. A general discussion of the laminar-flow problem 
aiong.with a summary of the triple-deck formulation and the resulting boundary- value 
problem is presented in the section "The Laminar-Flow Problem." The following sec- 
tion, entitled "The Numerical Method,” deals with all aspects of the numerical methods 
used to solve this problem. First, an inverse iterative scheme to solve the coupling 
between the various layers of the triple deck is discussed. Then finite- difference meth- 
ods based on the Keller "box" scheme (refs. 4 and 5) are formulated to solve the triple- 
deck equations. Methods are also discussed for the evaluation of the Hilbert integrals 
'arising from the analysis. Also presented in the section on numerical methods is a 
description of the accuracy and convergence properties of the numerical methods. A 
discussion of the results for both zero and nonzero incidence follows in the Results.^ Solu- 
tions to the symmetric problem have also been recently obtained by Jobe and Burggraf 
(ref. 6) and by Veldman and Van De Vooren (ref. 7). Detailed comparisons between the 
present solution and those of references 6 and 7 are provided. Also the solution for the 
drag of a finite flat plate at zero incidence is compared with the experimental data of 
Janour (ref. 8) and with finite- difference solutions of the complete Navier- Stokes equa- 
tions recently obtained by S. C. R. Dennis, who provided numerical data from his unpub- 
lished results; The solution for the velocity profile in the wake of the symmetric solution 
is compared with experimental data of Sato and Kuriki (ref. 9). 

Currently, understanding of turbulent interactions at trailing edges is rather less 
complete. Recent attempts to develop a rational theory of turbulent trailing- edge flows 
include the investigations of Spence (ref. 10) and Kvichman (ref. 11). In reference 10, the 
main correction to boundary-layer theory is assumed to arise from the pressure change 
across the wake generated by the singular curvature of the inviscid trailing streamline. 
This leads to a jet-flap model of trailing-edge flows. Spence's model is -inconsistent and 
leads to unacceptable oscillatory solutions downstream. The failure of Spence's theory 
is caused by the neglect of convective acceleration terms in the normal momentum equa- 
tions. It is interesting that a scaling analysis of Spence's interaction equation indicates 
the need to retain these terms in the lowest order theory. The investigation of Kiichmah 
in reference 11 is based on a Lighthill model of turbulent boundary layers near the trailing 
edge, which is treated as an inviscid rotational flow. Some examples of rotational flow in 
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a wedge-shaped compression corner are presented. Although these results are of inter- 
est, there is no attempt to develop a complete, rational theory and no consideration is ; 
given to the trailing-edge region on a lifting airfoil. . 

A more promising approach for a rational theory of turbulent flows follows from 
asymptotic expansions of the time-averaged Navier-Stokes equations in the limit of large 
Reynolds numbers. Similar asymptotic techniques have been applied to noninteracting 
turbulent boundary layers by Mellor (ref. 12), Yajnik (ref. 13), and Bush and Fendell 
(refs. 14 and 15) and to transonic shock-wave— boundary- layer interactions by Melnik 
and Grossman (ref. 16) and by Adamson and Feo (ref. 17). 

In the second part of the present investigation (Turbulent Trailing-Edge Flows) it is 
shown that asymptotic analysis leads to a three-layer description of turbulent interaction 
near trailing edges with a streamwise length scale that is on the order of a boundary- layer 
thickness. The flow in the outermost layer is governed by inviscid, linearized rotational- 
flow equations. The description near the wall requires two layers, involving just Reynolds 
stresses, in the middle layer and both Reynolds and laminar stresses in the innermost wall 
layer. The solution in the outer layer is unaffected to lowest order by the two inner layers 
and can, therefore, be completely determined independently of the details of the inner lay- 
ers. This leads to a Lighthill model for the outer problem that must be solved to deter- . . 
mine the pressure distribution and lift forces. Here, only the incompressible problem 
will be considered and a brief description of the essential features of the interaction model, 
together with a formulation of a boundary- value problem governing the outer inviscid flow 
will be provided. Also a sample solution for the skin friction determined from matching 
the inner and outer solutions is given and the results are compared with the low- speed 
data of Schubauer and Klebanoff (ref. 18). 

SYMBOLS • ; . 

parameters related to lift coefficient and constant defined by equation (19b) 

. displacement function in triple- deck theory 

constants related to behavior of A near the origin 

constant in asymptotic solution for |x| — <» 

constant defining y-grid distribution 

constant related to second- order boundary-layer solution for displacement 
thickness (turbulent) 


A 

Ai,Ajf 

Bp 
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D 

f 


Go 

hi,h2,h3 


H 


I 

Ip 


skin-friction coefficient, 




skin-friction coefficient of noninteracting boundary- layer solution at trailing 
e<^e (turbulent) 

constant appearing in asymptotic solutions of triple-deck equations, where 
i and j = 1, 2, 3, . . . 

lift coefficient 

coefficient of singular pressure gradient in wake 

constant related to singularity of inviscid solution near trailing. edge 

normalized shear stress gradient in computation plane, dr/dy 

constant in asymptotic solution for |x|— <» 

constant in triple-deck solution for drag coefficient. . 

shear stress gradient in Z- direction, 8 t/sZ 

Blasiu's function 

similarity function related to symmetric triple- deck solution for |x| — -<» 

Hakkinen-Rott similarity function related to triple- deck solution for X - 0 

functions appearing in differential equations for Hj, H 2 , and Hg, 
respectively 

shape factor, 6^62 

similarity functions related to triple-deck solutions for X — '->» 
integer, running index for X-mesh 
integer related to X-mesh 
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integers defining X- mesh 


integer, running index for y-mesh 


Finner number of points in inner mesh in wake 


number of mesh points in y-mesh on plate 


number of mesh points in y-mesh in wake 


parameter defining mesh distribution, where i = 1,2, 


je(y+) 


mixing length, turbulent flow 


L length of plate or chord, dimensional 

Lq,Lj,Ljjj 3 x parameters related to X- mesh 
p pressure 

p free- stream pressure, dimensional 

flO 

P normalized pressure appearing in triple-deck theory 

Pj,P 2 »P 3 constants related to behavior of P near origin 

R Reynolds number, UoqL/i^ ' 

Rjj functions related to asymptotic behavior of U as Z — «> , where 

i and j = 1 , 2 

Tj integrals defined by equations (45), where i = 1, 2, 3 

u velocity in streamwise direction 

u^- friction velocity, 

U normalized velocity in streamwise direction in triple-deck theory 


free- stream velocity, dimensional 
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velocity in direction normal to plate 


normalized velocity in direction normal to plate in triple-deck theory 
Coles wake function (turbulent flow) 

Cartesian coordinate along plate 

I 

normalized coordinate along plate in triple- deck theory 
end point of boundary- layer calculation 

computational coordinate normal to plate for laminar flow and physical 
coordinate for turbulent flow 

coordinate defining outer boundaries of y-mesh 

physical coordinate normal to plate for laminar flow 

normalized angle of attack in triple- deck theory 

angle of attack in radians 

constant appearing in behavior of triple- deck solution near origin 
, for a ^ 0 

Gamma function 

boundary- layer thickness 

displacement thickness 

momentum thickness 

boundary-layer camber, 

small parameter, for laminar flow and ^ for turbulent flow 

small parameter related to wall- layer thickness, (c^r) in turbulent flow 



5 vorticity (turbulent flow); also variable defined by equations (38) and (39) 

77 independent variable in similarity solutions of triple- deck equations, 

|Zl/2|2X|^/^^ 

e parameter used to scale wake location, related to wake centerline 

. K Karman constant, approximately equal to 0.41 

\ constant, equal to 0.33206, appearing in Blasius solution, f*'(0) 

Xj normalized skin friction at trailing edge in triple-deck solution 

AioAo 2>^23 integrals defined by equations (46) 

j ■ • ■ ■ . 

V coefficient of kinematic viscosity, dimensional 

^ variable used in definition of X-grid distribution 

IT Coles wake parameter 

p density, dimensional 

a parameter used to scale wake location, related to wake thickness 

T normalized skin friction in laminar study; also Reynolds stress in turbulent 

problem 

skin friction 

tf/ normalized stream function in computational plane 

stream function 

O) relaxation parameter 

Superscripts: 

* dimensional quantity 
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+ 


denotes wall-layer variable in turbulent analysis i 


^ denotes blending- layer variable in turbulent analysis 

' denotes differentiation with respect to indicated variable and perturbation 

quantity in turbulent flow 

Subscripts: 

B bottom surface of airfoil 

BL perturbation quantity arising from upstream boundary layer in turbulent 

analysis 

e local quantity evaluated at edge of boundary layer 

inv perturbation quantity arising from outer inviscid flow in turbulent analysis 

T top surface of airfoil 

TE quantity evaluated at trailing edge, X = 0 

THE LAMINAR- FLOW PROBLEM 

Problems of laminar flow at large Reynolds numbers are usually analyzed by a 
combination of inviscid- flow and boundary-layer techniques. This approach is based on 
asymptotic expansions of the complete Navier-Stokes equations in the limit of Reynolds 
numbers approaching infinity. The inviscid and boundary-layer equations arise as the 
basic equations governing the leading approximation in the outer and inner regions, 
respectively. This approach leads to accurate and useful solutions of viscous-flow prob- 
lems in many instances and has tended to dominate the history of fluid mechanics. How- 
ever, in spite of its central role in fluid floWs, the underlying structure of the asymptotic 
expansions are only relatively well understood for flows that are not separated and for 
geometries that are smooth. 

It is well known that the inviscid and boundary-layer descriptions break down near 
separation points or near singular points of the geometry, such as sharp leading edges, 
corners, and trailing edges. A comprehensive review of these matters, including a 
discussion of higher order approximations, has been given by Van Etyke (ref. 19) and 
Stewartson (ref. 20). 
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In the trailing-edge problem, the nonuniformity of the basic expansions, is caused 
jointly by a discontinuity an the surface boundary, conditions at the trailing edge and by a 
singularity in the inviscid solution at the trailing edge of a lifting airfoil. The disconti- o 
nuity in boundary conditions leads to a singularity in the boundary- layer solution at the 
trailing edge that is described by Goldstein's near wake- solution (ref. 21). Goldstein's 
solution shows that the displacement surface develops a sharp corner with a .vertical tan-r. 
gent on the downstream side of the trailing edge, as illustrated in figure 1. Goldberg and 
Cheng (ref. 22) have examined the second-order inviscid solution over a finite-length flat- 
plate at zero incidence and have demonstrated that the inviscid flow over this displace-., j 
ment surface is singular, with the induced pressures approachii^ plus (minus) infinity 
on the. downstream (upstream) side of the trailing edge. 

The inviscid solution for subsonic flow near the trailing edge of a cusped airfoil at 
angle of attack exhibits a square- root singularity in the surface pressure distribution,vas - 
sketched in figure 1. The surface pressure is bounded but the pressure gradients are; 
unbounded as the trailing edge is approached. The singular pressure, gradients lead to 
singularities in the boundary-layer solution and to breakdown of the basic asymptotic -, 
expansions. , J 

There have been numerous attempts to correct these defects and to develop an 
asymptotic theory that is uniformly valid at trailing edges, but for the most part, these* ' 
were completely unsuccessful. It was only in the recent work of Stewartson (ref. 1) and,' 
independently, Messi ter (ref. 2) that a correct a:nd rational treatment of the flow near 
trailing edges was given. In these works it was shown that the flow develops a character- 
istic multilayer structure near trailing edges that also arises in many other laminar 
interaction problems and which is referred to by Stewartson as a triple- deck model. A 
general discussion of viscous problems involving triple- deck structure is found in the ‘ 
recent review by Stewartson in reference 20. . . 

The Triple- Deck Formulation : - . 

Stewartson and Messiter presented a rational treatment for the flow near a trailing 
edge. By using the method of matched asymptotic expansions, these investigators have 
shown that solutions of the Navier- Stokes equations near the trailing edge can be developed 
in asymptotic series in the limit of large. Reynolds numbers. The solutions are cast in 
terms of a fundamental small parameter e given in terms of the Reynolds number 

e = R"l/® ■ . , (1) 

where R - is the Reynolds number based on the length of the plate and the constant flow 
velocity far from the plate. Stewartson and Messiter considered the idealized case of 
incompressible flow over a finite flat plate at zero incidence. The theory was extended 
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to angle of attack by Brown and Stewartson in reference 3 for angles of attack a* on the 
order Of e^/?. Although. the theory of references 1 to 3 was developed for Incompressi-' ; 
ble flow over finite flat plates, the basic theory is applicable to more general airfoils 
provided the airfoil is closely approximated by a flat plate near the trailing edge. For 
example, this requires the trailing-edge angle of an airfoil with thickness to be less 
than o(e2). The equations for the leading approximation are also applicable to coihpress- 
ible flows provided the Mach number of the inviscid solution is not near one in the trailing- 
edge region. In these cases compressibility effects enter only through a scale transfor- 
mation Of dependent and independent variables as given in references 3 and 20. 

The structure of the triple-deck region is sketched in figure 2. The flow upstream 
and downstream of the triple-deck region is governed by standard inviscid and boundary- 
layer equations. The leading term in the outer inviscid region is given by constant uni- 
form flow, while the solution in the upstream boundary layer is given by the Blasius solu- 
tion and in the downstream wake by a modified Goldstein near wake solution, as described 
in reference 3. In the intermediate region between the Blasius and Goldstein region, the 
flow develops the multilayer structure sketched in figure 2. The ratio of the length scale 
of each region to the length of the plate L :is also indicated in the figure. The stream- 
wise length scale is o(e^), which is an order of magnitude larger than a boundary-layer 
thickness. Viscous effects are important only in the lower deck where the solution is 
governed by classical (incompressible) boundary- layer equations. Both pressure and 
viscous forces are negligible in the main deck to lowest order. The main role of the 
essentially passive main deck is to transmit flow deflections generated by the sublayer to 
the outer edge of the boimdary layer. These flow deflections provide an inner boundary 
condition for the solution in the upper deck which is governed by inviscid small- disturbance 
equations. The solution in the inviscid upper deck is governed by elliptic partial differ- 
ential equations which provide for the long upstream influence that was missed in many 
previous theories. 

' From the preceding discussion, it can be. seen that the triple-deck formulation leads 
to a description of the flow as an interaction between the outer inviscid stream and the 
displacement thickness generated by the sublayer. The solution in the inviscid upper 
deck can be reduced to an integral relationship between the surface pressure and the flow 
deflection generated by the sublayer. 

Solution of the triple-deck problem is thus reduced to that of determining solutions 
to the boundary- layer equations valid in the sublayer. These solutions must match the 
rotational flow in the main deck and must result in a displacement thickness and pressure 
distribution that satisfies the linear integral relationship arising from the outer solution. 

The notation employed in references 1 to 3 and 20 varies. Here the notation 
employed by Stewartson in reference 1 will be followed with some exceptions. Physical 
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quantities are denoted by an asterisk, free- stream quantities by the subscript and the . 
plate length by L. The quantities x* and y* are Cartesian coordinates parallel and 
transverse, respectively, to the plate with origin at the brailing edge, u* and v* are 
velocity components in the x*- and y*-directions, respectively, p* is the pressure, 
p the density, a* the angle of attack, and \ is a constant, equal to 0.33206, associated 
with the Blasius solution for the wall shear stress. Nondimensional variables for the 
lower deck are given by 


X = \®/^x*/e3L Z = \^/^yyeSL 

U = u*/e\l/'^U«, . V = v*/e3\3/4u^ 

P = (P* - 




( 2 ) 


(These scalings are for incompressible flow; for compressible flow, see refs. 3 v . 
and 20.) The Reynolds number R is given by where u is the kinematic vis- 

cosity coefficient. 

For future reference, the solution for incompressible, inviscid flow over a flat plate 
of length L at incidence ct* is given by (for y* = 0) 


V* = 0 


u* = 


X* + B 


[(-x*)_(L + X*)] 


1/2 


sgn y" 


(-L <x* <0) 


(3a) 


u* = - V* = 


.* : + B 


x*(L + x*)^^ 


1/2 


(0<x*).: (3b) 


The lift coefficient Cl corresponding to the solution is given by 


Cl = 2ira*(l - ^ 


(4) 


where B ' is a constant to be determined. The value of the constant B, determined ' 
by the Kutta condition applied to the trailing edge, would be zero. Here, however, the 
trailing-edge interaction leads to a nonzero value which gives a viscous correction to the 
lift coefficient. In reference 3 it was shown that this constant is ole^L) which: leads to 
a viscous correction to the lift that is 0(e3). A nondimensional circulation constant aj 
is introduced according to the definition 


ai = 


_ B\ 


5/4 


e^L' 


(5) 
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The boundary- layer equations governing the flow in the lower deck are written in 
the form 


M + 9V 

ax az 


= 0 


( 6 ) 


ax az dX 92;2 


(7) 


where the pressure P is a function of X alone. 

Equations (6) and (7) are to be solved subject to the following boundary conditions: 


U - Izl + . . . 

(X - -«■) 

(8 a) 

U = V = 0 

(Z = 0; X < 0) 

(8b) 

P^(x) = Pg(X) 

IIV 

s 

(8c) 

fA.j.(X) 

(Z - +CO) 

(8d) 

U lz| + . . . — o/ 



1:Ab(x) 

(Z--CO) 

(8e) 


*v • • . ... 

where P-j.(X) and PgCX) are the pressures on the top and bottom of the X-axis, 
respectively, and A>p(X) and Ag(X) are perturbations of the displacement thickness 
from the undisturbed Blasius value at the trailing edge. 

Finally, the pressure distribution must satisfy the following asymptotic condition in 
. order to match the upstream inviscid solution 

Prji^g — -o;j/-X sgn Z (X — -®°) (9) 

where the subscript T,B is introduced for convenience to represent either the top of 
the bottoni and sgn Z should be taken as plus for Z > 0 and minus for Z < 0. The 
pressure should decay to zero for Z — +«> in order to match the Goldstein solution 
downstream . 

- Equation (8a) is a requirement that the velocity profile match to the inner portion 
of the Blasius solution far upstream. Equation (8c) follows from the requirement that the 
pressure be continuous across the wake. The application of this condition at the trailing 
edge (X = 0) is equivalent to the Kutta condition and serves to determine the constant aj 
defined by equation (5) and, therefore, the viscous correction to the lift coefficient. 
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Solution of the inviscid flow equations in the upper deck leads to the following inte- 
gral relationship between the functions Aip^gCX) and Fij- gCX): 


dAjCX) 


dAgCX) 

dX 



(Ida) 


(10b) 


where the double slash on the integral sign means that Hadamard "finite part" of the 
divergent integral is to be taken. The integrals in equations (10) do not converge in the. . 
ordinary, sense because of the behavior of the pressure for large negative X indicated 
in equation (9). , A form more suitable for computation is given in the numerical methods 
section. Equations (10) are an inversion of the relationships given in reference 3 and 
are in a form that is most suitable for the numerical procedures used in the present 
investigation. 

The preceding formulation indicates that to solve the triple-deck problem, the 
boundary- layer equations must be integrated subject to the vortical outer boundary con- 
ditions given by equations (8d) and (8e). The vorticity arises from the boundary -layer 
solution valid in the upstream flow. The upstream vorticity leads to additional algebra- 
ically growing terms in equations (8d) and (8e), as will be discussed later. 

The pressure and displacement functions appearing in the boundary- layer equations, 
and boundary conditions are unknown and must be determined as part of the solution of 
the boundary-layer equations such that the linear relation given in equations (10) is sat- 
isfied. The form of the pressure for large negative X is given by 

^T,B “ sgn Z (11) 

where aj is related to the lift coefficient by equations (4) and (5). The lift coefficient 
can be obtained by solving the boundary-layer equations and extracting the constant aj 
from the expansion given in equation (11). , 


Asymptotic Properties of the Triple Deck 

Equations (4) to (11) provide a complete formulation of the triple-deck problem. 
Useful asymptotic results were provided in references 1 to 3 and are extended and sum- 
marized in the following. (Corrections to a number of the signs are incorporated;, see 
also ref. 6.) 
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For X — the solutipn must approach the perturbed Blasius solution and have 
the form 


U = |z| + 1.784lxr^F^(7j) +,a|xl^/®Hi(7j)sgn Z + a2H^(T?) + a3|xr^/®H^{Tj) sgn Z + . . . 


and 



3l2X| 



(12a) 

(12b) 


where primes denote differentiation with respect to 77 . The similarity functions Fj^^^(tj); ' 
Hj^(fj), H 2 ( 7 J), and 113 ( 77 ) satisfy ordinary differential equations given in references 1 
and' 3. The differential equations, boundary conditions, and asymptotic behavior required 
in the present work are listed in the appendix. 

Other results for X — are 

P(X) = -0.34331x1"^^^ - 0.6867bilxr®^^ - 0.0816X"2 ln|X| + + . . . 

|x|2 

- o(|x|^/^ - ajlxl"^/^ +‘0.2368|xr^/® + . . .)sgn Z + . . .. ' ' (13^) 


A(X) = 0.32651x1“^ + aCi 2 lx|^^® sgn Z - ln|x| + In xj + a2(C2i + C 22 ) 

+ 0f3c33|x| sgn Z + ... 


(13b) 


M 

dZ 


z=o 


= 1 + 0.3106|x|”'^^^ - 0(2.1539)1x1'^^® sgn Z + 


(13c) 


where K and the Gjj's are known constants listed in the appendix. 

For X — +00 the velocity profile approaches the inner solution of Goldstein’s near 
wake solution. The behavior of P(X), Aq>^B(X), and U(X,0) are given by 


P(X) = 0.1717X'2/3 + 0.3433biX"®/?,- 0.0816X‘^ In X + d^X"^ + . . 


(14a) 


= 1.416(ix)^'^^ sgn Z + 1.416Q)^'^^ b^X"^/^ sgn Z - o 


I X^/^ + 2a^X^/2 




+ . . 


(14b) 
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(14c) 


u(x,o) = i.eiix^/^ + i.eiibjX"^/^ + o.oszx"^ + . . . 


The solution of the triple- deck equations develops a singularity at the trailing edge 
that is described by the Goldstein solution (ref. 21) for &= 0 and by the Hakkinen-Rott 
near wake solution (ref. 23) for a 0. In both cases the velocity profile has the follow- 
ing form for X and Z — 0 



where n is given by equation (12b) and Gq satisfies an ordinary differential equation 
studied in reference 23 and listed in the appendix. 

For zero incidence, the local solution for X — 0 is given as follows; 

For X < 0, 

P = P^E + PtE^ + (P 2 ln|X) + Pg)X^ + . . . (16a) 


A = + Ai^jjX - (s^/ysJColxl®'^^ + 1 AjX2 + . . . 


(16b) 


au 

az 


z=o 




(16c) 


and for X > 0, 


P = P^E + I CoX^/^ + P^X + . . . (17a) 

A = Arjig + AipE^ + ^3 ^^^/iO^CqX®^^ + ^ A^X^ + . . . - (17b) 

U = 1.40610^^^/^ + . . . ■ ' (17c) 

Co = 0.4089\|/^ (17d) 


The local solution for X — 0 when o; #0 is slightly more conaplicated for X < 0. 
Equation (16a) now takes the form 
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and 


dU 

dZ 


Z=0'^pc=0 


- ^1,T -az 


Z=0 ^=0 


(P3 )t,bx]x -2 ln|x| + . . . 

(18a) 

d Ag 



(18b) 


(18c) 

= ^l.B 

(18d) 


For X > 0, equation (17a) has the same form. Equation (17b) becomes 

At = (Ate)t + + ( 33/2/10)00X^/3 + i(Aj‘)Tx3 + . . . (18e) 

» -^B “ (^Te)b (^Te)b^ " + 2(^1 + . . . (18f) 


and equations (17c) and (17d) assume the forms, respectively, 

U = ycJ/Sjl/S ^ 

Co = Co(^1,tAi,b) (18h) 

where y and Cq are now functions of ^i^t ^1,B a-re determined by solving 
the merging asymmetrical shear flow problems of Hakkinen and Hott. 

The two sets of constants ^a^,b^„d|^ and (PjE’ ^TE» ^2» ^3» ^TE» ^E» 

A^, A^ , and g appearing in the preceding expansions are not determined by the 
local solution. The first set' relates to the for- field solution while the second set relates 
to the loc^ solution near the trailing e(^e. Values for all but d^ have been estimated 
for a; = 0 by fitting the asymptotic forms to the numerical solutions obtained in the 
present study. 

For future reference the behavior of the velocity profile for large Z and fixed X 
is given as 
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U = |z| + sgn Z + Agi ln|z| + ln|z|)sgn Z 

+ At,b sgn Z + ((ffCii/2Kl/2)AT^B + («y2Kl^^)|cn(C2l - C22) 


14 sgn ^ 1 


1/2 


1-3 


, , ,+ KC34 sgn Z> jZl - a2|Zj + 


(19a), 


The constants Cy and K are listed in the appendix and 


a2 = 1.784 


r(4/3) 

325/6(..y, 


(19b) 


The preceding expansion was obtained by expanding the solution of the boundary- 
layer equations for large Z and matching to the asymptotic expansion of the initial pro- 
file defined by equations (12). (See results in appendix.) This matching enabled the set of 
arbitrary constants appearing in equations (10) to be identified with the constants given 
in the appendix. Equations (19) are used to set the far- field boundary conditions in 
a finite difference procedure described in the following section. 


THE NUMERICAL METHOD 

. The boundary- value problem to be solved for the trailing- edge solution is illustrated 
in figure 3. The boundary-layer equations must be solved such that the solution matches 
the Blasius solution far upstream and the Goldstein solution far downstream. No- slip . 
conditions must be satisfied on both sides of the flat plate and asymptotic boundary condi- 
tions must be satisfied on each side of the boundary layer and wake for jz| — «>. The 
pressure gradient appearing in the momentum equation and the displacement functions 
Aij.(X) and Ag(X) appearing in the outer boundary condition must be determined such 
that they satisfy the linear integral relationship imposed by the outer inviscid solution. 

In addition, the condition that the pressure decays to zero as X — +«> must be imposed 
in order to match the’ Goldstein solution. In the present approach, a fixed point iteration 
between the inviscid and boundary- layer equations is employed. The principal difficulties 
in the numerical solution of the boundary-layer equations are due to a singularity at the 
trailing edge and to a slow algebraic decay of the solution for [xj and |z| - «. These \ 
problems are treated by using asymptotic solutions to set the far-field boundary conditions 
at finite distances and to describe the singular solution near the trailing edge. A highly . 
nonuiiiform mesh distribution is also employed to obtain propeir resolution near the trailing 
edge and to allow for the slow decay of the solution in the far field. 
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■ The bovindary-layer equations are solved by the Keller-Cebeci (refs. 4 and 5) finite- 
difference scheme for parabolic partial differential equations. This method is well suited 
to the present problem since it is second-order: accurate, unconditionally stable, and per- . 
mits highly nonuniform mesh distributions. The Hilbert transformations in equations (10), 
which provide the inviscid solution, are evaluated by special quadrature formula on the 
same mesh distribution employed in the boundary- layer calculation. This avoids the need 
to interpolate between the inviscid and viscous solutions. The resulting computation has 
uniform second- order accuracy, including the far field and the singular point at the trail- 
ing edge. 


The Iteration Scheme 

The iteration procedure is indicated in figure 3. The path of the iteration is in a 
direction inverse to that usually employed in similar viscous interaction problems. Here, 
the displacement functions A-j. g(X) are obtained from solutions of the inviscid equations 
(i.e., eqs. (10)) with the pressure distribution prescribed. The pressure distributions 
Pip g(X) are determined from solutions of the boundary- layer equations. Since the 
unlmown pressure, gradient appears in these equations, an additional relation is needed to 
complete their solution. This is supplied by the previous evaluation of the displacement 
function Afp^g(X) which provides an outer boundary condition for the solution of the vis- 
cous equations. This indirect iteration sequence is followed because it provides a con- 
venient and simple treatment of the trailing-edge singularity. In a conventional iteration, 
the solution of the boundary-layer equations for a prescribed pressure distribution results 
in a discontinuity in the slope of the displacement function at the trailing edge. This, in 
turn, leads to unbounded pressures in the inviscid solution and to divergence of the itera- 
tion sequence. , ■ . 

The iteration starts with estimated pressure distributions Pp(X) and Pg(X) 
which appear in the integrands of the Hilbert integrals. The integrals are evaluated by 
a second- order- accurate quadrature scheme to yield expressions for dAip^ciX and 
dAg/dX. The displacement functions Arj, and Ag are then obtained by integration 
using a trapezoidal rule with initial values determined from the upstream asymptotic 
expansions given in equation (13b).' This half- cycle yields an intermediate solution for 
the displacement functions A>j'^g(X). 

In the next half- cycle the boundary-layer equations are integrated. A minor diffi- 
culty arises because of the presence of the unknown pressure gradient in the differential 
equations. To deal with this problena the momentum, equation is differentiated with respect 
to ,Z. . This eliminates the pressure but increases the order of the equations from a third- 
to a.fourth-order system of partial differential equations. An additional boundary condi- , 
tion is required to close the system. This is supplied by using the known functions 
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Af b(X) in the asymptotic expansion given in equations (19) to yield a condition on the 
streamwise velocity component as |z| — «>. This condition and the conditions that the 
shear stress approach one for |z| — «> and that both velocity components vanish on the 
plate result in a well- posed problem. A finite- difference scheme, described in the next 
subsection, is employed to integrate the boundary-layer equations starting from an initial 
station far upstream of the trailing edge. Profiles at the initial station are determined 
from the first five terms of the asymptotic solution given in equations (12). The boundary- 
layer equations are solved by marching downstream to the trailing edge on the top and 
bottom of the platO independently. A local solution describing the singular behavior at 
the trailing edge is obtained by numerically solving the similarity equations, first consid- 
ered by Hakkinen and Rott in reference 23. The similarity solution is used to construct 
a "composite" profile across the sublayer at a station just downstream of the trailing 
ec^e. The solution is then marched downstream, employing two boundary conditions on 
each side of the wake, as indicated in figure 3. After completion of the sweep, the pres-' 
sure gradient is determined from the momentum equation evaluated on the X-^is at 
Z = 0. The pressure is then computed from a trapezoidal integration of the gradient 
Two arbitrary constants of integration, the trailing-edge pressure and the circulation 
constant aj (see eq. (11)) are evaluated by matching the pressure to the upstream data 
and by requiring the pressure difference to vanish at the trailing edge. 

.The. boundary- layer solution cannot be continued downstream to very large distances 
because of the -appearance of a growing solution P = PqX^/^ for X — +<». The solution 
is. induced by a wake thickness distribution + Ag) « AqX^/^ which appears in toe 

outer boundary conditions. The constant Pq vanishes and toe unwanted solution is 
excluded if the constant Aq is exactly equal to the Goldstein value (Aq = 0.892,. . .). 
However, because of toe finite arithmetic carried in toe computer, this solution cannot be 
excluded from toe numerical solution and it eventually dominates toe far-field behavior. 
Consequently, the boundary-layer solution must be terminated at a station X = Xp that 
is taken to be upstream of the region where toe spurious growing solution starts to dom- 
inate. This raises a minor problem, since a solution for the pressure distribution P(X) 
along the entire X-axis must be supplied for toe evaluation of the Hilbert integral. This 
is easily remedied by using an analytic expression to represent the pressure distribution 
downstream of toe terminal point X = Xp. In' the computer program an expression is ■ 
employed that matches both toe pressure and pressure gradient at X = Xp and has toe 
correct asymptotic behavior for X — +«>. Numerical experiments, to be discussed at 
the end of the section, have indicated that this procedure provides a smooth continuation 
of the solution downstream of X = Xp and has a negligible effect on the upstream 
solution. 

This half- cycle results in a complete solution for the pressure distribution which 
can be substituted into the Hilbert integral to obtain new estimates for the displacement 
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functions. The' integration is continued until the solution converges to a required toler- 
ance. As usually required in this type of problem, the solution must be underrelaxed in 
order to obtain convergence. 

The value of P(X) is relaxed according to the formula 

. P(X) = U.P(X)„^„ + (1 - o.)P(X)„,^ . 

where PqJjj is the pressure at the start of the boundary-layer computation, Pnew 
the pressure computed at the most recent sweep, and to is a relaxation parameter (to < 1) 
that is adjusted to obtain convergence. In the present scheme it is found that the value 
of to must be reduced as the extent of the streamwise interval is increased in the down- 
stream direction. Accordingly, the following strategy is employed. The calculation 
starts with a given relaxation parameter small enough to obtain convergence with an ini- 
tial choice of Xp = 3. The solution is converged with these choices, the terminal point 
is moved further downstream, co is reduced, and the calculations repeated. Converged 
results have been obtained starting with values of u = 0.15 and Xp = 3 and ending 
with 0 ) = 0.02 and Xp = 20.791. 

Solution of the Boundary- Layer Equations 

The boimdary- layer equations are solved by the Keller-Cebeci ’Tsox" scheme (refs. 4 
and 5)i The unknown pressure gradient is eliminated from the boundary-layer equations 
by application of a Z- derivative to the momentum equations. A stream function ’J' is - '- 
introduced and the boundary-layer equations are written as a system of four first-order 
partial differential equations. The wake thickness and centerline position become 
unbounded as X — «>. To control the wake growth the Z- coordinate is scaled such that 
wake position is bounded in the computational plane. A scale transformation is defined 
in terms of two parameters ct(X) and 0(X) by the following relation: 


where y 
a(X) and 


y = 


cKX) 


( 20 ) 


is a scaled coordinate in the direction normal to the plate and the functions 
6(X) are given in terms of A>j>^b(X) by 


<t(X) = 


o ■ -^b(^i)] 

0(X) = - -i|^Lp(X) + Ag(X^ + 


X < L,) 

(21a) 

All 

X 

(21b) 

X 

IIV 

(22) 
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where is a constant scale parameter and Lj is a small positive number that identi- 
fies the streamwise station where the wake solution is initialized. The functions oi(X) 
and 0(X) control the wake thickness and position, respectively. The choice of 0(X)' ' 
was previously discussed in reference 3. This scaling minimizes the variations of wake 
location in the far field and as a result, the computations can be carried out with fixed 
outer boundaries in the y-plane. In addition the coordinates approach similarity variables 
appropriate to the Goldstein solution for X — +«>. Scaled dependent variables are also 


introduced according to the relations 

a{X)^UX,y) (23a) 

U = a(X)u(X,y) (23b) 

T = t(X,y) (23c) 

D = a(X)~^d(X,y) (23d) 


where U, t, and D are, respectively, the stream function, streamwise velocity 
component, shear stress, and derivative of shear stress with respect to Z ^i.e., ^ ~ 
and where >p, u, t, and d are, respectively, scaled versions of the stream function, 
streamwise velocity component, shear stress, and derivative of shear stress with respect 
to Z. With these transformations and with the elimination of the pressure gradient, the 
governing equations can be written in the form 


aii/ 

— I- = u 
ay 


(24a) 


(24b) 


■^= d 

ay 


(24c) 


f =a(X)3(u^'-d^)-2a(X)2a'(X),/<i 


(24d) 


In this formulation the boundary condition on V given in equation (8b) is replaced by the 
equivalent condition on the stream function 

0 (y = 0;X<0) (25) 

A general, nonuniform rectangular mesh is introduced and equations (24) are dif- 
ferenced according to the box scheme along the lines indicated in figure 4. The X- columns 
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are labeled by an index I and the horizontal rows by an index J starting from J = 1 
on the lower boundary and continuing to J = Jp on the upper boundary, where Jp is 
the number of points in the y-mesh on the top of the plate. The first three equations do 
not involve X-derivatives. These are central differenced about a midpoint of a y-intervai 
on the most forward marching column. Equation (24d), which is nonlinear and involves 
derivatives in both directions, is central differenced about the midpoint of the box, as 
indicated in figure 4. The nonlinear coefficients are evaluated as four point averages at 
the midpoint of the box. This difference approximation leads to a nonlinear set of differ- 
ence equations for the vector unknowns («^,u,T,d)j along the column I + 1. The differ- 
ence approximation is second-order accurate and implicit since it couples all the unknowns 
along the I + 1 column. 

The boundary conditions along the plate involve the specification of the two compo- 
nents X < 0. CXiter boundary conditions are imposed on the 

vector components u and d. The conditions are given in the form of a ratio at the outer 
two points of the mesh Jp and Jp_i. The ratios for the conditions on top of the plate 
are computed from the asymptotic far-field expansion given in equations (19) as follows: 


o(X)ujp - E„| 

(Z ,X) 
( ‘’p J 

r 

Ap(X) Ri2(Zjp.X) 


^ *'P-1 


1 - Ag,(X) Rj^2| 

f2j 

1 “lp-1 } 




where 


(26a) 


(26b) 


Rjj ( ZpC ) + Z + Z^/2 + In Z ‘+ k ^/ 2 z " 1/2 in Z 


K 


1/2 


Rl2(ZP^) = 


aC 

|2K 


^ Ay(X) + - C 22 ) + Z 


1/2 - a2Z-3 


(27a) 


(27b) 


R2i(zp0’ = Z"3/2 + a2c^^z-2 . | in Z 


(28a) 


R22(ZP^) - ^^3^ ^2 ®CiiAt(X) + 0|3^-Cj^C22 + KC34 - | C1JC21) 


Z‘3/2 - 12a2Z‘3 (28b) 
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where Z is related to the computational coordinate y by . . 

Z = fl(X) + a(X)y ' (29) 

This procedure for satisfying the far-field boundary conditions was motivated by the 
work of Ackerberg and Phillips (ref. 24). Similar expressions are applied to the bottom 
boundary of the mesh. For the wake computations (x 5 Lj), only the outer boundary 
conditions are to be satisfied and the expressions are identical to those of equations (26), 
(27), and (28). 

The difference equations are solved by a Newton- Raphson technique. The nonlinear 
equations are linearized about a previous estimate to form a linear system of algebraic 
equations for the perturbation quantities. ( 6 i^, 6 u, 6 t, 6 d)j . The differential equations 
and boundary conditions result in a linear system that has a block tridiagonal form. In 
the present problem the main blocks are. 4 x 4 square matrices. The equations are solved 
by. an efficient Gaussian elimination technique as described in reference 5. The form of 
the outer boundary conditions given in equations (26) automatically falls into this block 
structure. 

Solutions at the most recently computed station (e.g., station I in fig. 4) are employed 
as initial estimates. Quadratic convergence was observed to occur with these starting 
values. The iteration was continued until a convergence criterion based on the relative 
error was satisfied at all mesh points. The criterion 



is employed, where f stands for any one of the dependent variables. 

The calculation proceeds by marching in the X- direction starting from an initial 
station X = Lq. Initial profiles are determined from the asymptotic solution given in 
equations (12). The similarity functions Fj(t)), Hj(? 7 ), H 2 (t?), and 113 ( 77 ) appearing in 
the asymptotic solutions are determined from a numerical integration of the two-point 
boimdary- value problems formulated in the appendix. These solutions are obtained with 
the same subroutine employed in the marching calculation. 

When the trailing edge (X = 0) is reached, a composite solution is formed to describe 
the initial wake profile a short distance (aX = Lj) downstream of the trailing edge. The 
composite profile is obtained from a coordinate expansion for X and Z — 0. It is written 
as the sum of an outer and inner solution less the "common part." The structure of the 
local solution is similar to Goldstein's near wake solution except for the presence of a 
singular, self-induced pressure gradient in the similarity, equation. Solutions were first 
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obtained by Hakkinen and Rott in reference 23 and are further discussed in references 1 
to 3 and 20. The wake initial profile is given by 


1/3 

U(Li,z) = Ut(0,Z) + 1(1 Lij Go(t?) - tZ 

(z > 0) 

(31) 

U(Li,Z) = Ub( 0,Z) + i| 

1/3 . 

Li) G'o(n) - Xi,bZ 

(Z <0) 

(32) 


where Ut and Ug are the velocity profiles at the trailing edge on the top and bottom 
of the plate, respectively, Gq is the similarity function describing the "inner" Hakkinen- 
Rott solution, T] is a similarity variable defined by equation (12b), and and ^ 

are the skin- friction coefficients at the trailing edge. (See eq. (18d).) Profiles for the 
variables i//, t, and d can be formed in a similar manner. The function Gq(t]) is 
determined as part of the solution by integrating the two- point boundary- value problem 
formulated in the appendix. The solution in the wake is then continued downstream, start- 
ing from the initial wake station (X = Lj) and terminating at a station (X = Xp) chosen at 
the start of the calculations, as discussed in the beginning of this section. 


The pressure distribution is determined after completion of the forward sweep by 
integrating the streamwise momentum equation 


^ d(X,0) 

dX ■ <r(X) 


+ a(X)2 


t(X,0) 


a^^(x,o) 

9X 


u(X,0) 


3u(X,0) u^(X.O) do^(X) 

8X ■ 2 dX 


(33) 


. All quantities on the right side of equation (33) are known from the most recent sweep. 
Separate equations hold on the top and bottom of the plate. The pressure is determined 
by a simple trapezoidal integration 


P.J. g(i + 1) - Pt,b(i) = [x(i + 1) - X(I)| 


\ ® )ui/2 


(34) 


where the pressure gradient is evaluated by averaging equation (33) oyer the stations 
X(I + 1) and X(I). Two constants of integration are required to complete the solutions 
given in equation (34). These are determined from the asymptotic solutions given in equa- 
tion (13a), which give 




_ 0.3433 0. 2368a 

5/6 


PB(X=Lj=a^Lo - 


o'^l _ 0.3433 ^ 0.2368a 


l/R 


2/3 .1 .5/6 


(35a) 


(35b) 
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Note that these relations involve the unknown circulation constant aj^ which must be 
determined before the solution can be completed. The circulation constant aj is deter- 
mined from a generalized Kutta condition as follows. By definition 


Pt(0) = Pt(X = l„) + 
Pb(0) = Pb(X = Lo) * 


dX I 
dX 


(36) 


If equations (35) are substituted into equations (36) and the Kutta condition that the pres-:. ■ 
sure is continuous at the trailing edge ^i.e., ■.Prj.(O) = Pb(0)) is imposed, the following... 
relations aire obtained: , , 



(37a) 


(37b) 


Equations (37) should be interpreted as asymptotic relations valid for Lp - The 
present results indicate that the solution is not overly sensitive to the magnitude of Lq 
and that z-i can be evaluated to two decimal places for Lq = -17. Equation (37a) has 
been employed in conjunction with equation (34) to determine the surface pressure by 
sweeping equation (34) from the trailing edge. 

The differential equations are differenced in Cartesian coordinates on a nonuniform 
mesh. The grid-point distribution is determined from simple transformations that map 
a uniform grid to a nonuniform grid.. The parameters of the .mapping are.ac^usted to. cpn- 
centrate: mesh points in regions of large gradients that occur on the axis (y = Q) and at the 
trailing edge (X = 0). The distribution of y-grid points is given as follows.for the upper 
half-plane; -’ The mesh points in the lower half-plane are obtained by refle ction. about the , 


X-axis. • .. . v’ : • ■ > ; , 

On the plate side (X < 0) the y(J) distribution is defined by the relations. ' 


y(J) = T 




M 


1 + Bp(l - ?) 


(38a) 
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and 




(38b) 


(IJJSJp) (38c) 

where is the upper boundary of the computational domain, J is a running index, 
and is a parameter employed to control the relative spacing of the increments. 
Equations (38) reduce to a uniform mesh for Kj = 1 and to a nommiform mesh with a 
concentration of mesh points near the wall for K| < 1. 

On the wake side of the field, a two-piece grid consisting of a fine uniform grid was 
employed near the axis, and a stretched grid was used in the outer part of the wake. The 
grid for X > 0 is defined by the following relations. For the inner region 

= Wr P - ^ •'inner! 0 ly t (39a) 

and for the outer region 

I'W* ' W * Pm ■ W) 1 . Bp(i - 0 [' * "3" ■ «“]' . <**'’> 

where 

^ ^Jinnw (Jinner < J = Jwl yinner .< y ^ yj^) 

where is the upper boundary of uniform mesh region, J is a running index, 

Jinner number of mesh points in the inner region, is the total number of mesh 

points employed in the wake, and K2 ^d K3 are parameters that control the mesh dis- 
tribution in the outer re^on. They are chosen such that (1) the mesh increment is contin- 
uous across the boundary between the inner and outer regions and (2) the outermost incre- 
ment y(j\v) ** y(Jw ' equal to the corresponding increment on the plate. The mesh 

distribution in tte marching or X- direction is chosen to provide for a concentration of 
mesh points on the wake side of the trailing edge and for a gradual stretching in the far 
field for X — ±<». 
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On the plate side (X < 0) the following transformation is enaployed: 


X(I) = Lo4 





(40a) 


where 


It> - I 

l = (l sisip; Lo sx so) (40b) 

Lq Is the value of X at the initial station, I is a running Index, ip is the number of'^ 
X-mesh points used on the plate, and K4 is a parameter that controls the relative spacing 
of the mesh points on the plate. Equations (40) generate a uniform mesh for K4 = 0 and 
a noniiniform ndesh for K4 > 0. A relation between K4 and the minimum mesh incre- 
ment Aj is given by ■ ' 

(Ip - 1) + (Ip - 1 )[ai(Ip - i)/Lof 

" ~ 1/3 

(ip- »)- [*i(ip-i)/U 

With equations (40) and (41), a mesh distribution can be generated with a specified mini- 
mum increment Ai at the trailing edge that smoothly expands to the initial point (x = Lq, 
I s'l). ' ■■ ^ ^ ; -u. 

The streamwise mesh distribution in the wake is given in three parts: a nonuniform 
region near the trailing edge that concentrates points near the. origin, an intermediate 
region with a uniform mesh, and a nonuniform region with an expwding grid in the down- 
. stream direction. ' A good distribution is, generated by the following relations: 


X(I) = 


(I, - Ip)* r (I, - Ip - 1)* 


X(I) = X(I - I) + 42’ 


ii:a), = x(i2) 


^2(1-13) 


i-K5(l-l3) 


3 


■; (lp <l gll) ’ ■ '‘ (42a) 

(^1 < I =^ 2 ),' ^2b) 
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(43) 


r" 

^2(13 + 1- 12) 

■jl/3" 

y 1 \ 

^'max. " ^(^2) 


\J3 + ^ " hj 

V 


J 



where A£ is a parameter that controls the first mesh increment. in the wake,-. I 2 , ■■ 
and I 3 are the values of the running, index I that separate the three mesh regions, 
and LjQax coordinate of the downstream boundary of the mesh. Note that I 3 is 

also equal to the total number of points used in the streamwise direction. Attention is 
called to the fact that L ^av “Ot equal to the terminal point of, the boundary-layer 

calculation, X = Xp. 


Evaluation of the Hilbert Integral 


The most recent sweep of the boundary-layer equations provides updated solutions 
for the pressure distributions required for the evaluation of the Hilbert, integrals given 
in equations (10). The main difficulties in the numerical evaluation. of these integrals are 
associated with the infinite range of integration, the algebraic singularity in the integrand 
for X the pole singularity at X = X<, and the infinite pressure gradient for 
X-0+. 

The first two problems are treated by dividing the integration interval into a nxunber 
of segments. The outer two segments contain the unboimded interv^ds X - +■» and 
X - -<». In these regions, the pressure distribution is approximated by the asymptotic 
expressions given in equations (13a), (13b), and (14a), and the integrals are evaluated in 
closed form. This reduces the numerical problem to one involving an inte^ation over 
a finite range and also provides for a correct evaluation of the singular "finite part" ' 
integral for X — The remaining inte^als are over a finite range and are evaluated ' 
by numerical quadrature using the mesh distribution employed in the boundary-layer cal- 
culation. Difficulties with the pole singularity are avoided by evaluating the integrals 
only at the midpoints of the ihte^ation intervals used in the quadrature. Excessive trun- 
cation error due to the sin^ar pressure gradients hear the trailing edge is avoided by 
using a special quadrature formula" that accounts for this behavior. Oh the plate, the inte- 
gration interval is split into two segments" -«« < X 5 Lq and Lq < X < 0, while in the 
wake, 0 < X < Lj, Lj < X < L 2 and L 2 < X < +«>. Note that segment boundaries Lj, 
L£, and L 3 need not line up with the boundaries of the mesh defined by equations (42). 
With this division, the integrals can be expressed as the following summation 

At,b(X) = f-Tj - sgn Z^Tg + T 3 + Ajq + Apg + A 23 ) (44) 
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where the far-field expressions have been used to evaluate the integrands of the inte- 
grals T2t and Tg as follows: 




o \/-XidXi pLp 
X-Xj “^ij. 


exi 


(-1/3)1 fLo 


£ 


dX, 


|/^(x - Xi) /3(6)2/3 (_x)5/6(x - Xj) 


T = - MM 

2 3\/3 *^-00 


dX^ 


(-x)*/“(x - Xj) 

T - 0.892 p*” 

3 S (x,f /3(X - X,) 

AloW - I" 






(45a) 

(4Sb) 

(4Sc) 

(46a) 

(46b) 

(46c) 


The integrals in equations (45) are evaluated in closed form and those in eqxiations (46) 
are evaluated by numerical quadrature. 

The range of integrations of the integrals in equations (46) is segmented by using 
the mesh distribution X(I) employed in the boimdary-layer calculation, and the integrals 
in equations (46a) and (46c) are expressed as a finite sum of integrals over the mesh incre- 
ments X(1 +1) - X(l). The individual integrals over these increments are then evaluated 
in closed form by using a piecewise linear approximation for the pressure distribution 
Pt q(X) and/or P(X) over the mesh increment. The integral in equation (46b) is 
evaluated with a piecewise linear approximation for the function P(X) - ^ c2/3x2/3, 
where Cq is the constant appearing in the Hakkinen-Rott similarity solution (eqs. (18)). 
With this procedure all integrations are second- order accurate. The displacement func- 
tions can be evaluated from a trapezoidal integration as follows; 


At,b(I + 1) = At,b(I) + |?(I + 1) - X(I^At3(i + 


(47) 
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with evaluated from equations (44) to (46).. Two constants of integration 

appear. These are evaluated from the asymptotic solutions given in equation (13b) at 
the initial station X(l) = Lq, which yields 

I 1-1 I il/6, o 

A'p — 1 ) “ 0*3265 |Lq| + 0£Oj^2|^o| Z ~ ck ^21 


ilnjLol + 103(21/3) 


+ a 


*(C 


21 


+ €22)+ “ c 


331 


- 1/6 


sgn Z 


(48) 


This procedure results in a convenient method for determining the values of Aij. g 
at the mesh points X(I) given the pressure at the same points. A number of nuinerical 
experiments have been carried out for sequences of mesh distributions and for pressure 
distributions that could be integrated in closed form. These results clearly indicated that 
the quadrature errors reduced quadratically with mesh size and that the preceding eval-. . 
nations of Aj b(X) were second-order accurate at all points of the mesh including 
those near the trailing edge. The results of this study will be presented in a separate 
publication. 


Accuracy and Convergence Considerations 

A number of numerical experiments were carried out to checlc the accuracy of the 
complete program. These tests were carried out for the symmetric problem (i.e., a; * 0) 
by using a version of the program that was modified to allow for the symmetry of the solu- 
tion. The angle-of-attack terms. were deleted and a symmetry condition was imposed on 
the wake axis. With the modified program it was necessary to cominite the solution only 
in the upper half -plane; thus the nimiber of mesh.points„ required was reduced by one-half. 

' ' . Calculations were performed to determine the effect of varying the loca:tions of the 

upstream (X = Lq) and downstream (x '= Xp) .boimdaries of the mesh and the position of 
the upper, boundary- yj^. The number of mesh points employed in the horizontal and verr 
tical directions were also varied, as were the parameters controlling the relative spacing- 
of the grid. Computations were carried out using up to 99 points normal to the plate and , 
300 points in the streamwise direction. These results indicated that. a good distribution 
of mesh points is generated with the following choices.- i. . ' 

For the y- mesh ' " . - • v 

= 8.0 Ki = 2/5 (49) 

and 25 points are employed on the plate side ^Jp = 25j and 44 points on the wake side 
= 44j of the trailing edge. Of the 44 points in the wake, 12 (Jinner “ 
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tributed uniformly in the region 0 s y 5 0.5555. With these choices the minimum incre- 
ment occurs on the axis and is equal to 

(^y)min = 0.1333 (X < 0) (50a) 

(^y)min = 0.0427 (X > 0) (50b) 

The mesh spacing smoothly increases to a maximum at the upper boundary where it is 
equal to (for all X) 

(^y)max “ 0.8027 (51) 

For the X-mesh the upstream Lq and downstream Ljjjax boundaries are chosen 
as 

Lo = -17 l^max = 30 (52) 

Fifty-one points are employed on the plate and 124 points in the wake. The boundaries of 
the various segments in the wake are taken at 

Ij = 81 I3 = 175 (53) 

and, therefore, a total of 175 grid points are employed in the streamwise direction. 

The initial station in the wake is taken at 1=57 which occurs at 

X(57) = 0.004136 (54) 

The minimum mesh increment on the plate occurs at the trailing edge and is equal to 

(^)min = 0.05 (X < 0) (55) 

In the wake the minimum increment occurs at the wake initial station X(57) and is equal 
to 

(AX)min = 0.002432 (X > 0) (56) 

These choices lead to very high concentrations of mesh points in the initial parts of the 
wake. About 10 points of the y-grid fall in the inner region of the wake initial profile 
where the solution is described by the Hakkinen-Rott similarity solution. Numerical 
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experiments have indicated that solutions on this mesh are accurate to about one part in 
the third decimal place. 

Convergence criteria were set on both iterative loops employed in the program to 
achieve three-place accuracy. The error tolerance used in the Newton- Raphson solu- 
tion of the difference equation (see eq. (30)) was set at 10“^. This resulted in solutions 

O 

to the difference equations that are accurate to 10' or better in two cycles per stream- 
wise step at most stations. A third cycle is occasionally required near the trailing edge. 

The overall interaction between the boundary- layer and inviscid programs is con- 
tinued until the solution has converged to the third decimal place. Twelve iterations were 
required to converge the main loop with w = 0.15 and Xp « 3. Each cycle consists of 
a sweep through the boundary layer and the evaluation of the Hilbert integral. Most of the 
computer time is taken in the boundary-layer routine. The computations were performed 
on an IBM 370/165 digital computer and required about 20 seconds per cycle, or about 
4 minutes to complete the first 12 iterations. A total of 28 additional cycles were 
employed to move the terminal point Xp downstream to Xp = 20.791. The influence of 
the value of Xp on the upstream solution was investigated and was found to amount to an 
increment of no more than 0.002 in the solution at X = Xp, which decreased rapidly 
upstream. Similar conclusions hold for the angle- of- attack problem except that the com- 
puter times were about doubled because of additional mesh points on the lower side of the 
flow field. ; 

The present program used significantly fewer iterations and less computer time 
than that required by similar methods developed by Jobe and Burggraf (ref. 6) to treat the 
zero angle- of-attack problem. This is apparently due to their, use of a fixed point itera- 
tion scheme to solve the finite- difference equations and to the use of a separate iterative 
scheme to compute the pressure at each streamwise station. Jobe and Burggraf were 
able to use larger values of w in the outer loop and, as a result, obtained solutions with 
somewhat fewer outer cycles. However, this advantage did not nearly overcome the 
longer cycle times required in their program. It should be pointed out, however, that the 
symmetric problem does not involve free parameters and, hence, needs to be solved just 
once. Therefore, computing efficiency is not a real issue in this problem. It was, how- 
ever, important to develop an efficient code for the full problem since there are twice the 
number of mesh points and since solutions must be obtained for various values of the nor- 
malized incidence a. 


RESULTS 

The triple-deck formulation reduces the trailing-edge problem to one involving a 
single parameter o', a normalized angle of attack defined in equations (2). The computer 
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program described in the previous section has been employed to obtain solution for two 
values of a equal to 0 and 0.10. Initial estimates to start the iteration were obtained 
from the approximate solution of Messiter (ref. 2) for a = 0 and the linear solution of 
Brown and Stewartson (ref. 3) for nonzero angles of attack. Numerical experiments car- 
ried out in the study indicates that the present solution for the symmetric case is accurate 
to three decimal places. The symmetric solution was obtained by using a special version 
of the code in which the angle- of- attack terms appearing in the boundary and initial condi- 
tions were set to zero. In addition, a symmetry condition was imposed on the axis (i.e., 
j//=T=0 at y = 0) and the solution was computed only in the upper half-plane. The ' 
solution for a = 0.10 was obtained with an early version of the code that employed a • 
somewhat coarser mesh and, hence, is likely accurate to just two decimal places. 

The results obtained for the symmetric problem are compared with solutions of the 
triple- deck equations recently obtained by Jobe and Burggraf (ref. 6) and Veldman and 
Van De Vooren (ref. 7). The computations in these studies were based on finite- difference 
techniques that differed in a number of respects from the method employed here. The 
main difference being that a second-order Crank-Nicplson scheme was employed in refer- 
ences 6 and 7, while a Keller-Cebeci box scheme was employed in the present study. The 
computations in reference 6 were carried out bn a nonuniform mesh using up to 180 points 
in the vertical direction and 480 points in the streamwise direction. In reference 7 a non- 
uniform mesh was employed with a maximum of 40 points in each direction. These are to 
be compared with the 24 x 175 point nonuniform mesh used in the present computation. 

The inviscid solution was obtained in reference 6 and in the present study from a numer- 
ical quadrature of the Hilbert integral. In reference 7 the inviscid solution was deter- 
mined from a finite- difference solution of Laplace's equation in the outer deck by using a 
40 X 40 point mesh distribution. The iterative techniques used in the present scheme 
appear to be more effective and require significantly less computer time than the methods 
employed in references 6 and 7. Comparisons given in this section indicate that the over- 
all agreement between the three sets of solutions is quite good, with differences amounting 
to a few parts in third decimal places at most points of the flow field. However, the solu- 
tions of reference 6 are somewhat less accurate near the trailing edge due to the poor 
resolution of the trailing-edge singularity obtained with the uniform mesh distribution used 
in that study. The computation in reference 7 loses some accuracy in the far field due to 
the large mesh spacing used in that region. 

Solutions to the symmetric problem are presented in figures 5 to 17. The pressure 
distribution on the axis is given in figures 5 to 7. The effect of the wake in generating a 
significant favorable pressure gradient on the plate is clearly shown in figure 5. The 
pressure starts from the free- stream levpl far upstream and falls to a value of 
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^TE “ 


(57) 

at the trailing edge. This should be compared with the values of P-j>g of -0.388 and 
-0.392 obtained in references 6 and 7, respectively. The pressure then rises steeply 
from this trailing-ec^e value to a small positive maximum and then approaches the free- 
stream value slowly from above. These results clearly show a large adverse pressure 
gradient in the wake just downstream of the trailing edge. The pressure gradient is 
bounded on the upstream side and is unbounded on the downstream side of the trailing 
edge. The numerical solution is seen to blend smoothly into the asymptotic far- field 
solution for X — and match smoothly with the singular solution at the trailing edge 
for X — O'*’. The coefficient of the third term of the trailing-edge solution (eq. (17a)) 
has been extracted from the present numerical results as 

i»l = -0.52 ^ (58) 

The agreement between the present solution and the solutions of references 6 and^ 7 is 
quite favorable, with the differences between the results being indiscernible on the scale 
of figure 5. 

The pressure distribution near the trailing edge is shown on a greatly expanded 
scale in figure 6. Differences between the three sets of results are apparent on this 
scale. The present results and those of Veldman and Van De Vooren are virtually identical 
with the three term expansion given in equation (17a) with the constant Pj given by equa- 
tion (58). The use of a nonuniform mesh with a fine grid near the trailing edge provides 
excellent resolution of the singular trailing-edge behavior in the present calculations and 
in those of Veldman and Van De Vooren. The results of Jobe and Burggraf, which were 
obtained on a uniform mesh, definitely appear to have a higher truncation error and to 
lose some resolution as the origin on the wake side of the trailing edge is approached. 
Their results, however, appear to improve as the mesh size is reduced. 

On the basis of analytical considerations, Stewartson in reference 1 has indica.ted . , 
that the pressure gradient is finite on the plate side of the trailing edge and that loga- 
rithmic terms must arise in the expansion of the pressure distribution as X — 0”. (See 
eq. (16a).) The present results plotted in figure 7 seem to confirm Stewartson 's conjec- 
ture. In figure 7, the numerical solution for the pressure gradient is compared with the 
analytic expression given in equation (16a). The numerical constants ^2’ ^3 

were extracted from the numerical solution and were found to have the following values: 

= -0.301 ?2 = 0.12 P3 = -0.14 (59) 
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This is to be compared with the value P^E “ "0.278 given in reference 6. The numer- 
ical solution clearly indicates a vertical tangent at the origin and shows creditable agree- 
ment with the analytical solution using the constants given previously. Four mesh points 
in the region X < -0.5 fall on the analytical curve given in figure 7. 

The solution for the skin friction is given in figure 8. The solution is seen to match 
smoothly to the weak interaction solution given in equation (13c). The results clearly show 
the strong effect of the wake-induced pressure gradient on the skin friction. The skin 
friction at the trailing edge is increased by a factor of Aj over the Blasius value, where 


Xi = 1.351 


(60) 


This is to be compared with the values of t^( 0) = 1.343 and r^(0) = 1.352 predicted 
in references 6 and 7, respectively. Comparisons of the triple-deck solution with the 
second-order boundary- layer solution of Schneider and Denny (ref. 25) are given in 
reference 6. 


The solution for the centerline velocity in the wake is given in figure 9, together with 
a comparison of the Goldstein solution for X — +« and with the Hakkinen-Rott solution • 
for X — 0". Both analytic solutions exhibit an X^/^ behavior and appear as linear dis- 
tributions in the scale used in the figure. Also included is a comparison with two terms 
of the Goldstein solution. The second term, involving the constant b^, corresponds to a 
: shift in the origin of the asymptotic solution. The value of the constant b^ has been 
extracted from &e present numerical solution, as will be discussed later in this section. 

It c^ be seen that the triple-deck solution provides a smooth blending between the trailing- 
edge and far-field solutions. The effect of the shift is clearly evident in the results. . 


In figure 10, the solution for the centerline velocity is compared with the results of 
references 6 and 7 on an expanded scale near the origin. The present results show good 
agreement with the solution of Veldman and Van De Vporen and with the singular solution 
of Hakkinen and Rott right to the trailing edge. The honuniform mesh ernployed here and 
in reference 7 permits a very high resolution of the singular trailing-edge behavior. The 
results of Jobe' and Burggraf (ref. 6) again show higherj truncation errors and somewhat ■ 
poorer resolution near the trailing edge owing to the larger mesh intervals employed in 
their uniform mesh solutions. 


The solution of Veldman and Van De Vooren employs a highly stretched mesh with 
relatively large mesh increments in the region away from the trailing edge. The results 
in figure 10 ihdicate that this leads to somewhat larger truncation errors in the down- 
stream region than those that arise in the present solution. 

The present solution for the displacement function A(X) is compared with the far- 
field asymptotic expansions in figure 11. Again the numerical solution blends smoothly 
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into the asymptotic solutions for X — The effect of the origin shift in the down- 
stream solution is again evident in the comparison. The inclusions of the bj term in 
the asymptotic solution is seen to extend the region of agreement well into the near field. 
The results are seen to be in good agreement with the solutions given in references 6 
and 7, with no apparent differences on the scale employed in figure 11. ; 

Comparison of the present results for the slope of the displacement surface with 
the solutions of Jobe and Burggraf shows some discrepancy as indicated in figure 121 
The dashed line in figure 12,',representii^ the solution of reference 6, was obtained from , 
a graphical reading of a figiure in reference 6 using an automatic digitizer. Some of the 
differences are surely due to errors in reading the graphical data. However, the main 
differences are in the trailing-edge region, and these are likely caused by the larger grid 
spacing used in reference 6. The present results clearly show the vertical tangent at the 
origin implied by the singular solution given in equation (17b). This is seen more clearly 
on the expanded scale used in figure 13. \ In figure 13,i the present numerical solution for 
A’(X) is compared with the singular expansion given in equation (17b) and with the solu- 
tion of reference 6 as tabulated in the thesis of Jobe (ref. 26). The present solution is 
seen to blend very smoothly with three terms of the singular solutions. The constants in 
equation (17b) were evaluated by fitting equation (17b) to the present numerical solutions 
and were found to have the following values: 

0.338 A!Ipg = 0.402 AJ = -1.3 Aj = -2.1 (6i) 

These are to be compared with the values Aijg = 0.335, A^g = 0.335, and Aj = 0.56 ’ 
given in reference 6. The results agree relatively well with the solution of reference 6 
except for the grid point nearest the trailing edge and the values of the constant A][. 

The constant bj^ appearing in the second term of the far-field solution has been 
evaluated by fitting equations (14b) and (14c) to the present numerical solutions for A(X) 
and U(X,0). The results are denoted by b^(X) and bu(X) and are displayed as func- ' 
tions of X in figure 14. The results seem to approach a limiting value that is given by 

bi = -0.285 ± 0.005 (62) 

which is to be compared with the value bj = -0.27 ± 0.03 quoted in reference 6. -Also 
shown for comparison is a similar plot taken from reference 26. The difference between 
the two sets of results is likely caused by somewhat higher truncation error and by the 
abrupt termination procedure employed in the calculation of references 6 and 26. 

The drag coefficient for the finite flat plate can be evaluated from an integration of 
the skin- friction distribution on the plate as follows: 
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Cq = 1.328R'1/2 + dgR"'^/® + o(r'^) (63a) 

where the constant ^2 is given by 

do = f [t(X, 0) - l]dX 

i 

Use of the present solution for the skin friction and 
ceding integral leads to the following evaluation: 

d£ = 2.660 ^ (63c) 

In figure 15 the drag coefficients predicted by equations (63) are compared with experi- 
mental data obtained in 1935 by Janour (ref. 8) for flow of oil over a finite flat plate. 

Also included is a comparison with solutions of the full Navier-Stokes equations recently ' 
obtained by Dennis in 1973 (as mentioned in the Introduction) for Reynolds numbers in the 
range 1 s R s 200. The results in the figure show that the correction to the Blasius 
resvilt is large in this range and that it is accurately predicted by equations (63) to within, 
a few percent for Reynolds numbers as low as R = 10. 

Dennis later extracted the value of d 2 by fitting an equation of the form of equa- 
tion (63a) to his numerical solutions. His results for d 2 are plotted in figure 16 together 
with the limiting values (i.e., for R — +<«) predicted by the triple-deck solutions as 
obtained in the present study and in references 6 and 7. The agreement of all three triple- 
deck solutions with Dennis* results is quite good with the present solution yielding the 
best agreement 

The preceding results indicate that the triple-deck solution is accurate over a sur- 
prisingly wide range of Reynolds numbers. Indeed the maximum difference with the 
Dennis solution for the drag coefficient is about 8 percent at a Reynolds number R of 1. 
The close agreement with the Navier-Stokes solutions implies that the next term in the 
asymptotic solution, which is formally on the order of o(r"^) must be very small. Fur- 
ther comparisons and discussions of the drag coefficient are given in references 6, 20, 
and 26. 

Sato and Kuriki (ref. 9) carried out wind-tunnel experiments on the flow inithe wake 
of a thin plate. The flow was determined to be two-dimensional. The plate was 300 milli- 
meters long and the flow velocity was 10 meters /second. The investigators were primarily 
interested in exploring the transition of the wake from laminar to turbulent flow. They 

value of d 2 = 2.644 attributed to the present 'authors in references 6 and 7 was 
obtained on a coarser mesh than the one employed in the present computations and is, con-r 
sequently, less accurate than the above value. 


(63b) 

a trapezoidal integration of the pre- 
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measured mean as well as fluctuating velocity profiles in the wake. The Reynolds num- . 
ber of the test was 2.1 x 10® (e = 0.216). The mean velocity profile given in figure 17 
was measured at a station 30 millimeters behind the trailing edge where the flow was 
fully laminar (the nonlinear transition region started about 40 to 60 millimeters behind 
the plate). At this Reynolds number the pressure peak predicted by the triple-deck theory 
occurs at 40 millimeters behind the plate. Thus, the measuring station for the profiles 
in figure 17 was in a region where the theory predicted a strong adverse pressure gradient. 

The triple-deck solution was used to construct a composite velocity profile at 
X = 30 millimeters. The solution was represented as the sum of an "outer and inner 
solution" minus the "common part," as follows: 

= f’(7?) + 0.4945A(X)|f"(77) - f"(0)| + 0.1642(Uinner * z) (64a) 

Uoo 

where 

Tj = 0.4941Z (64b) 

and the physical distance normal to the wake axis is given by 

ymin = 0-325Z (64c) 

The function f(tj) is the Blasius function for the semi- infinite flat plate solution and 
Uinner^Z) is the triple-deck solution for the wake profile at X = 2.49. The displace- 
ment fvmction at this station is given by 

A(2.49) = 1.052 (64d) 

The profile given by equations (64) is compared with the measured profile of Sato 
and Kuriki in figure 17. The theoretical and experimental profiles are seen to be in good 
agreement across the entire wake. The main differences occur in the outer region where 
viscous and pressure gradient terms have been neglected in the theoretical solutions. 

Also indicated is the centerline velocity predicted by the one- term Goldstein solution. 

The effect of the interaction in reducing the centerline velocity is significant and readily 
discernible in this experiment. It is also of some interest to call attention to the fact that 
transition was observed to start at a station 40 millimeters behind the plate, which coin- 
cided with the location of the theoretical pressure peak in the wake. This result suggests 
that self- induced pressures may play an important role in the transition of a wake from 
laminar to turbulent flow. As a corollary it also indicates that the effect of wake- induced 
pressure gradients may have to be accounted for in theoretical transition calculations. 
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The numerical methods developed in the present investigation have been found to 

provide an effective means for solving the triple-deck equations. Because of the overall 

efficiency of the differencing and iterative techniques employed, it is practical to use these 

methods to solve the angle- of-attack problem. An early version of the code described in. 

previous sections was used to obtain solutions for a (normalized angle of attack) equal . 

to 0.10. This solution was obtained without the terms of o(q! 2) or greater that appear : 

in the outer boundary condition given in equations (26) to (27) and in. the initial condition 

given by equations (12). These terms are quite small and are believed to have a small ' 

influence on the solution at the value a = 0.10 for which the computations are carried 
2 

out. 

In figure 18 the solution for the pressure distributions on the t< 5 ) and bottom of the 
plate and on the wake axis is compared with an approximate solution developed in refer- 
ence 3. The approximate solution was based on a simple linearization of the triple equa- 
tions about a linear streamwise velocity profile. As was noted in reference 3, the linear- 
ization is clearly not valid in the wake, where the velocity gradient 9U/^ must vanish 
on the axis. However, the errors in the wake are not expected to have a strong influence 
on the solution upstream of the trailing edge. If this holds true, the linearized solution 
should provide a reasonably good approximation to the angle- of- attack solution. The 
results in figure 18 bear this out. The agreement between the present numerical solu- 
tion of the full triple- deck equations and the linearized solution given in reference 3 is 
seen to be quite good. The effect of incidence on the pressure distribution in the wake is 
barely noticeable on the scale used in figure 18. The circulation constant aj appearing 
in the formula for the viscous correction to the lift coefficient in equations (4) and (5) is 
determined as part of the present solution and is given as 

aj = 0.55 (65) 

This is to be compared with the value determined from the linearized solution of Brown 
and Stewartson, namely 

ai = 0.79 (66) ^ 

The agreement for aj^ is not nearly as good as for the pressure distribution but is prob- 
ably as good as one should expect from such a simple approximation. 

In figure 19, the numerical solution for the pressure distribution is compared with 
the inviscid solution on the plate and with asymptotic solution valid for upstream. The 
numerical solution is seen to blend smootUy into the upstream asymptotic solution. Com- 
parison with the zero angle- of-attack solution given in figure 5 indicates that the approach 
to the far-field solution is much slower in the angle- of- attack case. The difference 

^Computations including all terms appearing in equations (26) to (29) have been 
carried out since this paper was written and the results confirm this conclusion. 
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between the numerical and inviscid solutions in the far field is due mainly to the circu- 
lation term aj. This term can be interpreted as a shift of the origin in the far-field 
asymptotic solution for X — -«>. The effect of the shift is clearly evident in the compar- 
ison in figure 19. The results in figure 19 clearly show the large changes in the inviscid 
pressure distribution induced by the wake. These changes take two main forms: one is 
the shift in the pressure distribution in the far field mentioned previously, and the other 
is the complete change in the shape of the pressure distribution near the trailing edge. 

The effect of the induced pressures on the solution for the skin friction is indicated 
in figure 20. In this figure, the numerical solution for the skin friction on the top and 
bottom of the plate is compared with the asymptotic solution valid far upstream. The 
numerical solution is seen to, join smoothly to the asymptotic solution. The far-field 
solution plotted in figure 20 contains three terms: the basic Blasius value, the second 
due to inviscid pressure gradients induced by incidence, and the third term due to the 
favorable pressure gradient induced by the wake interaction. These terms combine to 
yield a skin- friction variation that is quite small except very close to the trailing edge. 

The effect of incidence is seen to have a fairly large effect on the value of skin friction at 
the trailing edge, which was equal to Cf(0) = 1.349 for a = 0. Clearly, a much larger 
value of incidence will be required to drive the skin friction to zero on the upper service. 
The present result seems to indicate that the point of vanishing skin friction should first 
arise at a station upstream of the trailing edge. 

The solution for the displacement functions A-j<(X) and Ag(X) is given in fig- 
ure 21 where it is compared with the far-field asymptotic solutions. The strong effect 
of incidence in displacing the wake centerline is clearly evident in this result. Again, 
attention is called to the smooth blending of the numerical and asymptotic solutions in the 
far field. 

TURBULENT TRAILING- EDGE FLOWS 

The numerical solutions obtained in this study and in references 6 and 7 have com- 
pletely confirmed the triple- deck model and the local asymptotic solutions developed in 
references 1 to 3. These works provide a sound theoretical framework for analyzing 
laminar interactions at trailing edges of streamlined bodies. Unfortunately, the boundary 
layer on an airfoil usually undergoes transition to turbulent flow at the Reynolds numbers 
of interest. Most theoretical methods for predicting the effect of boundary layers on air- 
foil characteristics are based on classical second-order boundary- layer theory. Although 
not generally recognized, second-order boundary- layer theory for turbulent flows breaks 
down at airfoil trailing edges. Consequently, the theory does not provide a satisfactory 
basis for computing boundary- layer corrections to inviscid airfoil solutions. It is, there- 
fore, important to develop a systematic theory for treating turbulent interactions in airfoil 
problems. 
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In the present study, a formal asymptotic theory has been developed for turbulent 
interacting flows, following broadly along the lines of the laminar triple-deck theories. 

Of course, turbulent boundary layers differ greatly from laminar flows and are controlled ‘ 
by quite dissimilar physical mechanisms. Particularly significant is the fact that turbu- 
lent boundary layers tend to retain high velocities much closer to the wall than in a lam- 
inar’ flow. This effect is one of the main reasons for the small upstream influence 
observed in turbulent interactions. However, the general idea of the triple-deck approach ' 
in seeking formal asymptotic solutions of the Navier-Stokes equations in the limit of large 
Reynolds numbers is also useful for the turbulent, problem. 

Accordingly, solutions of the turbulent airfoil problem will be developed as formal 
asymptotic expansions of the full Navier-Stokes equations in the limit R — ■». The tur- 
bi:^nt analysis is carried out by using the basic framework developed in references 1 to 3 
for laminar trailing-edge problems. Incompressible flow over a thin airfoil with a cusped 
or nearly cusped trailing edge at angle of attack is considered. The airfoil is assumed 
sufficiently thin so that the inviscid flow is described to lowest order by the flat plate ' 
solution given, for example, by equations (3) and (4) of this paper. Here, however, the 
angle of attack a* is assumed to be 0(1) and not necessarily small. Transition is 
assumed to occur upstream of the trailing edge and the boundary layers are assumed to 
be fully developed turbulent flows in the trailing-edge region. More specifically, the 
velocity profiles in the noninteracting region are assumed to have a small defect form in 
the main part of the boundary layer and the usual logarithmic behavior near the wall. 

In this section, the main features of a turbulent interaction theory that can be devel- 
oped under these general conditions are outlined. First the behavior of boundary-layer 
theory near trailing edges is examined. It is shown that a singularity arises in the solu- 
tioh%hich causes’a-breakdbwn of the second-order theory near trailing edges. This anal- 
ysis will show that a Kutta condition for the second-order solution cannot be satisfied and 
that the lift correction cannot be determined. 

The failure of second- order theory is resolved by the introduction of a locat solution 
that correctly describes the flow near trailing edges. The local solution is shown to have 
a three-layer structure that superficially resembles the triple-deck structure of the lam- 
inar problem. However, in the turbulent problem the physical mechanisms leading to this 
structore and the basic equations holding in each region are very different from those 
arising in the laminar flows. 

In the present discussion, only a very brief description of the turbulent interaction 
theory is presented. The asymptotic structure of the local solutions governing the flow 
near the trailing edge will be described and the boundary-value problem that must be 
solved to complete the solution in the trailing-edge regiOn will be outlined. Also a simple 
solution for the skin friction in the trailing-edge region that follows from the theory is 
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presented. The skin-friction solution is compared with experimental data of Schubauer 
and Klebanoff (ref. 18) and with numerical solutions of the boundary-layer ^nations. 

Second- Order Solution 

The second- order solution for viscous flow over an airfoil is determined in the fol- , 
lowing steps: 

(1) Determine the inviscid solution for flow over the prescribed airfoil shape and 
compute the pressure distribution on the surface 

(2) Solve the boundary-layer equations with the pressure distribution obtained from 
the inviscid solution and compute the displacement thickness 

(3) Use the displacement surface to compute the equivalent source/sink distribution 
on the airfoil and wake surface and solve the inviscid equations with the source distribu- 
tion as a surface boundary condition. Alternatively, the displacement thickness can be 
used to form an equivalent airfoil shape that serves as a new geometry in the inviscid 
solution. 

This well-known procedure can be embedded in a formal asymptotic expansion for 
large Reynolds numbers. However, because of the presence of a singularity in the inviscid 
solution, this expansion is not uniformly valid at trailing ecfees. Consequently, the second- 
order solution cannot be completed and the boundary-layer corrections to the lift coeffi- 
cient cannot be determined. Previous airfoil calculations based on second-order boundary- 
layer concepts relied on numerical smearing of the trailing- edge singularity to obtain 
solutions. 

A singularity appears at the trailing edge in the inviscid solution for all lifting air- 
foils with a sharp trailing edge. For an airfoil with a cusped failing edge,. the pressure 
distribution on the surface exhibits a square- root behavior. This is illustrated in figure 22 
where the steps leading to the nonuniformity of the second-order theory are outlined. Near 
the trailing edge, the pressure distribution is given by 

P* = Pte ■ y ^ 

where L is the airfoil chord and is a constant that depends on the incidence a* 
and the shape of airfoil. For a flat plate (or a sufficiently thin airfoil) can be deter- 
mined from the solution given in equations (3), namely 

Ca=a* ^ (67b) 

and = p^ and U,j,jg = U^. 
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The behavior of the displacement thickness near the trailing edge can be determined 
from the momentum integral equation 


dSg + 6jdx*(H + 2) 


dUe 

U*dx* 



( 68 ) 


where is the momentum thickness, H is the shape factor, Ug is the streamwise .. 
velocity external to the boundary layer, and Cf is the skin-friction coefficient The 
displacement thickness 6^ is given by 



(69) 


With‘the external velocity evaluated from the inviscid solution (eqs. (67)), the second term 
in equation (68) is, unbounded at the trailing edge. Since the variation of skin. friction due . 
to pressure variations can be shown to be equal to the order of the pressure, change,. the. . 
singular acceleration term in equation (68) can only be balanced by the gradient of 
momentum thickness. Thus, equations (67) to (69) lead to the following expansion for. 
the displacement thickness as x* - 0; 



(x* < 0) (70a)' 

(x*<p) (70b) 


where tte subscripts T and B refer to the top and bottom of the airfoil, respectively. _ 
The. boundary .layer thickens on the, top and thins on the bottom of the airfoil due to the 
imposed pressure distribution. This leads to an equivalent camber distribution 6^ given 



■*/ ' . • 



' ' 

6g - " ®1,b) “ ®c(®) Bq;\/-x*/L .+ . • . ... 

' .(71a) 

where 

aS(0) = l[6;.T(0)-5;.B(0| 

‘ ■ ■ 

(71b) 


Bo = - i CojsJ .p(% * 

x*=0 

(71c) 

The slope of the equivalent camber distribution given in equations (71) is singular at the 
trailing edge. The antisymmetric part of the second-order outer solution is determined 


by computing the inviscid flow over this camber distribution. It follows from potential 
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flow considerations that the second-order solution for thie surface velocity must have the 
following behavior for x* — 0~: 



where and are the constants defined previously and Ar is an arbitrary 
constant. 

The first two terms in equation (72) arise from the inviscid solution and the last 
two are induced by the interaction. Notice that both interaction terms are singular. The 
logarithmic term is induced by the singularity in the camber distribution, while the Ar . 
term is an arbitrary homogeneous solution that satisfies Laplace's equations and the 
boimdary conditions. Ordinarily this term would be excluded by the Kutta condition, which 
requires the solution to be bounded at the trailing edge. It is clear from equation (72) that 
this condition cannot be satisfied for any value of the constant Ar. Thus, the second- 
order solution cahnot be completed and the boundary- layer correction to lift cannot be 
determined. It is curious that this conclusion, which follows from the simple analysis 
given previously in this paper, has gone unrecognized in previous viscous airfoil analyses. 

Interaction Theory 

The results given previously clearly demonstrate that the standard second-order 
boundary- layer theory is not imiformly valid at trailing edges. To develop complete solu- 
tions of viscous airfoil problems, the basic theory must be corrected to better account for 
the flow near the trailing edge. In the present investigations, the method of matched 
asymptotic expansions was used to develop formal solutions for the trailing- edge region. 
This approach is based on the time-averaged Navier-Stokes equations with a turbulence 
closure employii^ a turbulent kinetic equation and an algebraic length- scale relation. 
Solutions were developed in terms of a small parameter e which here is related to the 
friction velocity u* in the noninteracting region upstream of the trailing edge. That is. 



where Cf q is the skin-friction coefficient at the trailing edge, as determined from solu- 
tion of the noninteracting boxmdary-layer equations on the top of the airfoil. Asymptotic 
solutions are developed for R - <» or equivalently for e 0. The analysis follows very 
closely a similar theory developed for interactions between turbulent boundary layers and 
normal shock waves in reference 16. 
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In the brief discussion of the theory given here, only a general description of the 
main results of the analysis will be provided. In addition, the present discussion will be . 
limited to incompressible flows. A more complete discussion of the analysis leading to 
these results along with a simple extension to compressible trailing-edge flows will be 
provided in future publications. 

From the present analysis, it has been determined that the streamwise extent of 
the interaction region at the traUing edge is on the order of a boundary- layer thickness; 
that is. Ay* ~ 6* » eL. The flow near the trailing edge was found to develop a multilay- 
ered structure, as illustrated in figure 23. The solution upstream of the interaction is 
divided into standard potential flow and boundary-layer regions over a streamwise length 
scale 0(L). The solution in the boundary layer has a two- layer structure typical of tur- 
bulent flows; an outer wake-like region and an inner wall layer. The velocity profile in 
the outer region has a small defect form, which on the top of the plate can be written in 
the form 

-^ = 1 + ef(y*/6^pc*) (74) 

e ,T 

where Ug is the velocity at the edge of the boundary layer and . 6^ is the local 
boundary-layer thickness. The velocity profile in the inner layer is expressed in a law 
of the wall form 



(75a) 


y+ = y*u*^T/i' C75b) 

where u is the kinematic coefficient of viscosity. Similar expressions hold for the 
boundary- layer profiles on the lower surface. The solution for noninteracting turbulent 
boundary layers has been embedded in a formal asymptotic structure by Mellor (ref. 12), 
Yajnik (ref. 13), and Bush and Fendell (refs. 14 and 15). These authors have shown that 
the law of the wall and velocity defect profiles appear as the leading terms of an asymp- 
totic expansion for R — «>. Thus, fully developed turbulent boundary-layer flows can be 
viewed as limiting solutions valid in this limit. The present analysis should be considered 
as an extension of these works to the trailing-edge interaction problem. 

The solution in. the .interaction region develops the three -layer structure illustrated 
in figure 23. ,, 
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The three layers required in the solution are 

(1) An outer, basically inviscid rotational stream. This regipn includes most of the 
boundary layer and a part of the irrotational flow outside the boundary layer that is on^ 
the order of a boundary- layer thickness. 

(2) An inner wall layer that is a continuation of the wall layer from upstream. 

(3) An intermediate, or blending layer that occurs between the outer and wall layers^ 
The biending layer is thinner than the outer layer but thicker than the wall layer. 

The wall-layer thickness is defined in terms of the parameter I introduced by 
Mellor in reference 12 

Ay*=^eeL ■ (76a)^ 

■ i 0 / 

where 



In the outer layer both the Reynolds and viscous stresses are small compared to the 
inertia terms in the momentum equation. Vorticity is generated in the upstream boundary 
layers and is convected, unchanged, along streamlines in the trailing-edge region. This 
leads to a description of the flow as an inviscid rotational stream. A similar model was ’ 
first proposed by Lighthill in 1953 (ref. 28) for treating interactions of oblique shock waves 
with turbulent boundary layers at supersonic speeds. 

The flow in the inner layer is a local equilibrium flow in the sense of Townsend 
(ref. 29). To lowest order, the total stress (viscous plus Reynolds stresses) is constant 
across the layer and the solution is completely determined by the local skin friction. 

An intermediate region is required because of a mismatch that develops between the 
Reynolds stresses in the inner and outer layers. In the outer layer, the Reynolds stresses 
are. frozen at their upstream values. In the inner layer, the Reynolds stresses are in a 
local equilibrium determined by the wall friction because of thegS mall- scale structure of 
the turbulence in this region. As a result, a discontinuity in Reynolds stresses develops , 
between the inner and outer regions. This discontinuity is resolved by the blending layer. 
Solutions in the blending layer are governed by linearized boundary- layer equations that 
involve Reynolds stresses, but not viscous stresses. Turbulent closure models are 
required to complete the lowest order solutions in the wall and blending layer regions but , 
not in the outer region. Displacement effects generated by the two inner layers are small 
and do not affect the first few terms of the solution in the inviscid outer region. Thus, the 
leading terms of the outer solution can be determined without consideration of the flow 
near the wall. 
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Next the form of the expansions in each region is considered. The expansion param- 
eter e defined in equation (73) is equal to the friction velocity of the upstream flow on 
the top of the airfoil. It should not be confused with the previous definition of e used in, 
the laminar study (i.e., Cj^aminar “ 

Outer Layer 

The coordinate stretchings for the outer region are given by 

X = (x*/L)e-l y = (yVDe'^ (77) 

where the coordinate system employed in the laminar study is used. There will be some 
minor differences from the notation used for the laminar analysis but the changes will be 
clear and should not cause confusion. 

The solution in the outer layer is dominated by contributions from the irrotational 
airfoil solution. The first terms in the expansion are obtained by expressing the airfoil 
solution in the inner variables defined in equations (77) and expanding the result in powers 
of e. The next term comes from the upstream boundary layer and is determined by sim- 
ply adding the defect part of the upstream profile to the irrotational contribution. 

Physically, this approximation is based on the idea that velocity variation across 
the boundary layer is small in the limit of large Reynolds numbers. (See eq. (74).) Thus, 
the physical picture is one of a basically irrotational flow that is slightly perturbed by 
small shear flow disturbances near the wall. The resulting linear superposition of the 
boundary-layer and inviscid solutions leaves the pressure distribution in the field and the 
lift coefficient of the airfoil unchanged. The nonlinear interaction of these terms produces 
perturbations in the pressure distribution and lift coefficient. Thus, the expansion in the 
outer region is written in the form 

H 1 // = y + el/2j//.^^(X,y) + e)^gL(y) + + . . . (78a) 

Ute 

= u = 1 + el/2uij^y(X,y) + eugL(y) + e3/2u'(X,y) + . . . (78b) 

Ute 

where p and u are the nondimensional stream function and the streamwise velocity 
component, respectively, is the velocity at the trailing edge predicted by the irro- 

tational outer solution, and * denotes the corresponding dimensional quantities. The 
velocity components u and v are related to the stream function by the usual relations 

U = dp/dy V = -dp/dX (79) 
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The first two terms in equations (78) arise from the airfoil solution written in inner var- 
iables and then expanded in powers of e. For a general airfoil with a cusped trailing edge 
this yields the result ' 


’^inv = “ cos | 0 (80a) 

'^inv = Cojr^/^ sin i 0 (80b) 

where r and 0 are polar coordinates, with 0 measured from the positive X-axis, 
given by 

r = ^x2 + y2 (81a) 

0= tan-1 (y/X) (81b) 

Substitution of the assumed expansion into the time-averaged Navier-Stokes equations 
leads to the following equation for the perturbation stream function: 

v2,//’ = -?’(X,y) = - (82) 

dy2 

2 I 

where ^ is the Laplacian operator and ^ is the perturbation vorticity. This is a 
simple Poisson equation that relates the disturbance stream function to the perturbation 
.in vorticity The vorticity perturbation arises from the convection of the vorticity in 
the upstream boundary layer along the curved streamline of the irrotational airfoil solu- 
tion, as illustrated in figure 24. 

The Poisson equation must be solved subject to the boundary conditions 


xl/{X,y)^0 (r-<=o) (83a) 

v’ = - = 0 (y = 0; x s 0) (83b) 


where v' is the perturbation velocity normal to the surface. The perturbations in 
streamwise velocity and static pressure are given by the relations 



p' = -u' = 



(84a) 

(84b) 
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These boundary conditions lead to well-posed boundary- value problems for the outer 
solution. The outer solution leads to a "slip" velocity on the surface that is resolved by 
the inner layers. 

The solution of the boundary-value problem can be represented as the sum of a 
particular solution plus a complementary solution. The particular solution satisfies 
equation (82) but not the boundary conditions. This solution leads to a downwash on the 
surface that violates the boundary condition on v' (eq. (83a)). The complementary solu- 
tion is a solution of Laplace's equation that cancels this downwash. The solution for the 
particular integral depends on the form of the initial velocity profile UBL(y) and on the 
expression for the irrotational stream function general profiles, the partic- 

ular and homogeneous solutions must be found by numerical means. However, if the ini- 
tial profile is represented by a Coles law of the wallAaw of the wake correlation, a closed- 
form expression for the particular integral can be found by analytic function theory. Coles 
form for the defect profile on the upper surface can be written in the form 


= \ 


iln(y/6T) -^W(y/6T) 






(0 ^ y g 6^) 

' (^T<y) 

where 6,p is a nondimensional boundary-layer thickness defined by the relation 


(85) 


6^ — €L6rp 


(86) 


K US the, Karman. constant, - 5frp is the Coles wake parameter, and is the wake 

function which can be represented by a simple polynomial approximation to the cosine 
function usually employed in this description. Similar expressions hold for the profile on 
the lower surface. 

The particular integral evaluated in this fashion leads to closed-form expressions 
for the downwash velocity on the top and bottom of the plate. The homogeneous solution 
which cancels this downwash can be found in the usual way from thin airfoil theory. This 
leads to a representation of the homogeneous solution in terms of a Hilbert integral that 
must, in general, be evaluated by numerical quadrature. A Kutta condition, requiring the 
solution to be finite at the trailing edge, is imposed as part of the solution of the homo- ■ 
geneous problem. This condition determines the value of an arbitrary constant appearing 
in the trailing-edge solution that is directly , related to the circulation constant Ar 
appearing in the second-order boundary- layer solution (eq. (72)). Matching the trailing- 
edge solution to the second-order solution valid outside the trailing-edge region leads to 
an expression for the lift coefficient in the form 
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(87) 


Cl = 2ira*(l + In e + a 2 e^ + . . . j 

where aj is a known constant and a 2 is a constant to be determined from a com- 
plete solution of the trailing-edge problem. The constant aj as determined by match- 
ing to the upstream solution (eq. (72)) is given by 



where 6 it,b ^T,B a^re the (nondimensional) displacement thicknesses and shape 
factors at the trailing edge, respectively, as determined from the noninteraction boun(tory- 
layer solution. The nondimensional displacement thicknesses Siq' b are defined by the 
relation 

^1T,B “ ^^^®1T,B (®®) 

where Si-p g are the dimensional displacement thicknesses. Notice that the leading 
correction to the lift coefficient, as given in equation (87) is o(e2 in e). This term is 
completely missed in standard second-order theory which leads to a correction that 
is o(e2). 

Equations (87) and (88) indicate that the lift correction is due primarily to the differ- 
ence in boundary-layer thicknesses on the top and bottom of the airfoil. In the usual situ- 
ation, rp > 6j g, and the effect of the log term is to reduce the lift coefficient. The 
effect of this term is most important on rear-loaded airfoils where the rear loading tends 
to dramatically increase the difference in boundary- layer thickness on the top and bottom 
of an airfoil. The effect of the shape of the boundary- layer profiles also influences the 
lift correction through the values of Hrj. and Hg appearing in equation (88). Rear load- 
ing tends to make the boundary layers less fvill on the top of the airfoil compared with the 
bottom. This implies that Hrp > Hg and this, in turn, also leads to a reduction in lift. 
This effect is formally of higher order since 


Hq-^g = 1 + 0(e) (R - «>) (90) 

However, in practice, H is significantly different from one at Reynolds number of 
interest (e.g., H = 1.4 for a flat plate at R = 10®) and this effect can be numerically 
significant. 


Inner Layers 

Next briefly the form of the expansions in the two inner layers near the wall is con- 
sidered. (See fig. 23.) Only the solution on the airfoil surface upstream of the trail ing 
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edge is considered. The evolution of the inner layers into the wake leads to a similar 
structure. However, solutions in the inner layers of the wake are much more complex 
and have not yet been fully developed. A major uncertainty in the wake solution is con- 
cerned with the choice of a closure hypothesis to properly deal with a change in the sign 
of the Reynolds stress near the axis. 

Only the solutions on the upper surface are dealt with explicitly. The expressions 
to be presented also hold on the lower side with an obvious change of notation. In the 
innermost layer the solutions for the streamwise velocity u* and Reynolds stress t* 
are represented in the form of a law of the wall as 

u* = u+^T^(X;e)/p(y+,X;e) ■ (91a) 

T* = T*(X;e)T+(y+pC;e) (91b) 

= (91c) 

where T^(X;e) is the skin friction and u"*", T+, and y"^ are nondimensional wall-layer 
variables. 

Substitution of these variables into the time-averaged Navier-Stokes equations with 
€ — 0 leads to the well-known condition that the total stress (laminar plus Reynolds 
stresses) is constant across the wall layer to all orders in e. This conclusion follows 
from the fact that the wall- layer thickness is transcendentally small in e. The wall-layer 
formulation is completed by the choice of a closure condition relating Reynolds stress to 
mean velocity. Analysis indicates that a balance of production and dissipation in the wall 
layer is a rational result that follows from the turbulent energy equation in the limit e — 0. 
With the usual model of dissipation this leads to a mixing length formula 

' T+ = £2(y+)(au+/ay+)2 (92) 

where £(y'*') is the mixing length distribution. In this formulation the choice of jf(y'^) 
is strictly empirical. Careful consideration of the magnitude of the pressure gradients 
in the present problem indicates that the choices commonly used for moderate-pressure- 
gradient flows are appropriate here. For example, the two-layer model of Cebeci and 
Smith with a Van Driest damping factor is known to give very accurate solutions in incom- 
pressible wall layers. It is known that the mixing length distribution is linear for large y+ 
and that this leads to the usual logarithmic velocity profile for y+ — <». Thus, if 

£(y+) = Ky+ (y+-.oo) (93) 
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where k is the Karman constant, it follows that 

u* - ^T^(X;e)/p(K"l In y+ + gf) (y+ - «) (94a) 

T* • T^(X;c) (y+ « oo) (94b) 

where is a universal constant, independent of the local value of skin friction. With 
T*(X;e) expanded in a series in e, the preceding result leads to an asymptotic solution 
for e — 0 and y*" — » that must be matched to the inner limit of the outer solution. 

The solution for the velocity in the outer inviscid region has a similar behavior for 
y — 0. However, because the Reynolds stresses are frozen in the outer region, the coef- 
ficient of the logarithmic term in the outer solution remains constant, equal to the friction 
velocity in the noninteracting boundary layer upstream. Thus the solution in the wall layer 
does not match the outer solution. It follows that an intermediate or blending layer must 
be inserted between the outer and wall layers in order to obtain a continuous solution for 
the Reynolds stresses across the boundary layer. 

The requirement that the Reynolds stress be continuous leads to the condition that 
the shear stress term must be retained as a leading term in the streamwise momentum 
equation in the limit 6-0. This condition determines the thickness of the blending layer 
to be o(e2L). Thus a new stretched variable y is introduced to represent the solution 
in the blending layer, where 

y* = e2Ly (95) 

Consideration of the form of the velocity profile in the upstream region and in the outer 
and wall layers leads to an assumption for a solution in the blending layer of the form 


u = 1 + el/2uinv(X,0) + e In e(lA) + 

+ e2/2 In eu2i(X) + e2/2u22(X,y) + . • ^. (96a) 

V = e3/2^2(X,y) + . . . (96b) 

P = P^E + e^/2p.^v(X) + e3/2p2(X) + . . . (96c) 

t = e2|l + el/2T2(X,y) + e3/2 ln€T3j(X) + + • • *j (96d) 
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where uj^y(X,0) and ii 3 L(y) are deduced from equations (80b) and (85) and are given 

(97a) 
.(97b) 

The In terms appearing in the preceding expansions are required in order to 
match similar logarithmic terms in the inner expansions of the outer solution. These 
terms enter the expansion from the logarithmic behavior of the initial profile for small y. 
Consideration of the transverse momentum equation leads to the conclusion that the pres- 
sure is constant across the blending layer to the order considered in equation (96c). Thus 
the unknown function P 2 (X) appearing in equation (96c) can be identified with the surface 
value of the pressure distribution of the outer solution. The requirement that the blending 
layer solution match the solution in the wall layer leads to the conclusion that the expan- 
sion for the Reynolds stress must have the same form as the expansion for toe streamwise 
velocity component. It also follows that the relative change in skin friction in the trailing- 
edge region is on toe order of toe pressure change; that is, = o(e^/2) a,s indicated 

in equation (96d). 

Substitution of the preceding expansions into toe momentum equation leads to the 
simple result 

U2i(X) = -K-luinv(X,0) (98) 

The, solution for toe vertical velocity component V 2 is obtained from integration of the 
continuity equation * . . . • 

, - V2^^,y) = -|dUinv(X,0)/d3^y . (?9) 

r The solutions for U 22 T 2 are governed by two coupled first-order, linear, 

partial differential equations. These equations are derived from toe streamwise momen- 
tum equation and, from a turbjd.ent plQSure. hypothesis relating Reynolds stress to mean 
velocity. Since toe pressure gradients on an airfoil can be relatively large, the closure 
assumption is based on toe turbulent energy equation. Since the solution of these equa- 
tions is not being considered in this presentation, they will not be written out here. Note 
simply that the momentum equation leads to a balance of linearized convective terms with 
shear stress and pressure gradient terms. Only the Reynolds stresses contribute to the 
shear stress gradients. The pressure gradient term is impressed from the outer invii^cid 
solution (dP 2 (X)/dx). 


Uinv(X,0) = sgn y 

= «“^jin(y/6rp) - 27?^ 
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The turbulent energy equation leads to a balance of advection, production, and dis- 
sipation of the Reynolds stress perturbation T 2 - The contribution from pressure diffu- 
sion and other terms in the energy equation are formally smaller than these main terms. 


The present formulation leads to a simple form for the skin friction that is inde- 
pendent of the particular closure model assumed in the analysis. Matching of the blending 
and wall-layer solutions for u* given in equations (94) and (96) leads to the following 
expression for the skin- friction' coefficient: 


Cf(X) 

^f,o 


- Cp(X;e) + 2(^lnf^Cp(X;e) + . . . 


( 100 ) 


where Cf ^ = 2€^ and Cp(X;e) is a pressure coefficient based on reference conditions 
from the inviscid solution at the trailing edge. Equation (100) is a simple relation for the 
skin-friction coefficient in terms of the pressure coefficient. It is a direct requirement 
of matching and follows simply from the three-layer structure of turbulent boundary 
layers near trailing edges. It involves only two turbulence parameters Cf^o 
The parameter Cf is the skin coefficient upstream of the trailing edge and k is the 
Karman constant which enters from the logarithmic term in the initial profile. 


Although the skin- friction result was derived here in the context of the trailing-edge 
problem, it can be given a more general and useful interpretation. The three-layer stuc- 
ture of turbulent flows also appears to apply to situations with large imposed pressure 
gradients. In this case equation (100) is valid with identified with the skin-friction 

coefficient upstream of the large-pressure-gradient region. The pressure coefficient is 
then defined with respect to reference quantities at the beginning of the pressure change. 
To check these concepts the skin- friction coefficient predicted by equation (100) was com- 
pared with data of Schubauer and Klebanoff (ref. 18). In reference 18 Schubauer and 
Klebanoff measured the skin-friction and pressure distributions in turbulent boundary 
layers approaching separation in a moderately large adverse pressure gradient. The 
skin-frictiori coefficient was computed by using the experimentally determined pressure 
coefficient in equation (100). The results are compared with the data in figure 25. Com- 
parisons with the turbulent boundary computations of Bradshaw, Ferris, and Atwell 

(ref. 30) obtained with a turbulent energy approach are also included. 

. . ■ t ' 

In this figure the combination 1 - Cp(X;e) is referred to as the first term. This 
one- term solution is equivalent to the assumption of a constant local skin- friction coeffi- 
cient (i.e., a skin- friction coefficient based on the local dynamic pressure at the edge of 
the boimdary layer). This result correctly indicates the main trend of the skin-friction 
variation with pressure but is in relatively poor agreement with the data. The inclusion 
of the logarithmic term in the solution greatly improves the agreement with the experi- 
mental data. 
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The numerical solutions given in figure:. 25 were obtained with and without a curva- 
ture correction. The present results should be compared with the numerical solution 
without the correction since it is not included’ ih'the present results. The agreement with 
the uncorrected numerical solution is seen to be quite good throughout the pressure rise. 
The better agreement with the experimental- data and with the corrected numeric^ solu- 
tion is probably fortuitous. However, the comparisons in figure 25 clearly indicate .that- 
the logarithmic interaction term is large and that very good results are obtained with its 
inclusion. These results also suggest that' equation (100) can be made the basis of a sep- 
aration criteria. The comparisons given tend to confirm the multilayer structure of tur- 
bulent flows proposed in the present study. , 


Concluding Remarks 


In the present investigation a formal asymptotic description of turbulent interactions 
at airfoil trailing edges was developed. The most important result of the present study 
was the formulation of a boundary-value problem that governs the solutions for the inter- 
action pressure distribution and lift coefficient on an airfoil. The interaction can be 
described as an inviscid rotational flow governed by a linearized Poisson equa:tion. Work 
is currently in progress to complete the solution of these equations. Also the analysis 
has recently been extended to compressible flow and it has been detnonstrated that the 
basic formulation applies to the case with minor modification provided the Mach number 
near the trailing edge is less than one. Both of these developments will be described in 
future publications on the subject. 


The results discussed in this section are concerned with the lifting, basically anti- 
. symmetric problem. The effect of wake- induced displacement pressures, which was so 
important in the laminar problem, is absent to the order of the small parameter e con- 
sidered so far in the turbulent problem. The present solutions can be carried to higher 
order in e. The next terms in the series are likely to the order of e2 in e and e2. 
These terms involve the thickness effects of^,the boundary layers on the airfoil and in the 
wake. However, the resulting problem is symmetric and so the solution would not affect 
the lift coefficient to this order. The major unsolved aspects of the problem concern the 
structure of & Oxpansiohs M th¥itmer layers of the' wake. Further analysis- is required 


to clarify the nature of the solution of this complex problem. 



' APPENDIX 

THE SIMILARITY FUNCTIONS 

‘ The function F^(?7) arising in the symmetric problem for X large and negative 
satisfies the differential equation and boundary conditions 

Fj” - 36(rjF[ - F^) = (Ala) 

Fj(0) = Fj(0) = 0 F'{iv)^0 (fj^cc) (Alb) 

It follows that 

Fi - [3-'^/6(1/3)I/(-1/S)]77 . („ ^ «) (Ale) 

The functions Hi(i7) arising in the angle- of- attack problem for X large and nega- 
tive satisfy the differential equations and boundary conditions 

Hi".- IStjZh-' + 9(4 - i)(7jH[ - Hi) = ^ (A2a) 

Hi(0) = H[(0) = 0 H[’(7])-0 (A2b) 

where 

hj = 9/2^/® ' (A2c) 

h2 = (3/2^/^)(3HiHi - nf) (A2d) 

h3 = (3/2^/^)(2 HiH 2 - hJHi + SHgHi) (A2e) 

It can be shown that the Hi(ti)'s have the following asymptotic behavior for 

Hi = |Ciitj3/2 + Ci277 + C13 + . . . (A3a) 

H 2 = C 21 T 7 In Tj + C 22 V + C23V^^^ + C24 +.. • . (A3b) 

H3 = C3 j7?^/2 In t? + C32 In tj + C33TJ + C3477I/2 + C35 + . . . (A3c) 
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APPENDIX - Concluded 


where the constants have been determined from analytical studies in references 1 
to 3 and from numerical solutions obtained in the present investigation to be 


Cii = -3(2)^/®(-2/3)!/(l/6)! 

“ “3*255 
Cm = C?j/2K 
C33 = -17.408 


Ci2= 6^/®(-l/3)! 
C22 ~ 3.082 ^ 

C31 = C„Cji/K 
C34 s 16.900 


Ci3 = -1/K 
C23 = C11C22/K 

C 32 = C 12 C 21 /K. 

C 35 = (C 2 I + C 22 )Ci 2 /k 


V • 


>’ (A4) 


J 


Where K = 3(2) 


1/3 


The function Gq(t}) governing the trailing-edge behavior satisfies the following 


differential equations and boundary conditions: 



' G”’ + 2GqG;; - G|? = 27(2)^/^Co 


(A5a) 

Gq(tj) — I 8 A .2 — Dfj" 0 

. (77 « 00) 

(A5b) 

Gq(tj) - I 8 X 1 gTj = Dg - 0 


(A5c) 

where j and X^^q are the values of the skin friction 

8U/az|^ and au/az|g 

at 


the trailing edge and Cq is a constant to be determined as part of the solution. Further 
details, together with typical solutions, are given in references 3 and 23. 
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Figure 1.- Trailing-edge flow. 
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Figure 2.- Triple -deck structure. 
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Figure 3.- Boundary -value problem for triple-deck solution. 



Figure 4.- Finite-difference approximation ("box" scheme). 
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Figure 5.- Surface pressure distribution; of = 0. 



Figure 6.- Pressure distribution near origin; a -0. 
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Figure 11.- Displacement function; . a = 0. 



Figure 12.- Displacement surface slope A'(X); a = 0 









Figure 15.- Drag of finite flat plate as a function of Reynolds number 

at zero angle of attack. 
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Figure 16,- Comparison between triple-deck and Navier-Stokes 
solutions for drag constant. 
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Figure 17.- Comparison of triple-deck velocity profile in wake of 
symmetric plate with data of Sato and Kuriki (ref. 9). 



X 

Figure 18.- Comparison of present solution with linearized solution 

of reference 3; o! =0.1. 
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Figure 21.- Displacement function; a. = 0.1. 
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Figure 22.- Second-order boimdary- layer theory - turbulent flow. 
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Figure 23.- Flow field structure - turbulent interaction. 
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Figure 24.- Outer problem. 




X(FEET) 

Figure 25.- Comparison of solution for skin friction with data 
of Schubauer and Klebanoff (ref. 18) and a turbulent 
boundary -layer solution. 1 ft ^ 0.3048 m. 
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